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Abstract 


A problem  is  considered  to  determine  the  tracking 
capabilities  of  an  estimator,  applied  in  optical  sensing. 

The  estimator  tracks  the  centroid  of  a one-dimensional 
Gaussian-shaped  intensity  based  on  time-space  point  process 
measurements.  The  centroid  assumes  to  move  dynamically 
as  a First  Order  Gauss-Markov  process.  Filter  performance 
is  described  by  steady-state  upper  and  lower  bounds  on 
mean-square-error  (MSE)  which  are  evaluated  as  a function 
of  two  physically  motivated  parameters:  average  number  of 
photons  detected  in  a coherence  time  of  centroid  dynamics 
and  mean  square  value  of  centroid  jitter  normalized  by  the 
square  of  the  beam  width.  The  parameters  establish  regions 
of  operations  where  upper  and  lower  bounds  converge  to  the 
actual  MSE.  Results  from  the  Monte  Carlo  simulation  demonstrates 
the  bounds' usefulness . Noise  measurements  from  dark  current 
or  background  radiation  are  included  in  the  simulation. 

Results  show  that  the  filter  is  very  sensitive  to  these 
measurements,  resulting  in  very  poor  tracking.  Ad  hoc 
methods  of  filter  tuning  and  residual  monitoring  are 
employed  to  improve  tracking  performance;  results  indicate 
that  filter  performance  can  be  improved  substantially 
through  residual  monitoring. 
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Fundamental  Limitations 
of  Optical  Trackers 


I.  Introduction 


Alignment  problems  exist  in  any  optical  system.  To  image 
a point  object  over  large  distances  requires  aligning  an  optical 
receiver  with  the  object.  When  the  object  is  within  the  receiver’s 
field-of-view,  an  image  is  formed  on  the  photodetector  surface. 
This  is  shown  in  Fig.  (1).  Because  the  optical  receiver  has  a 
relatively  narrow  field-of-view,  aligning  the  receiver  becomes 
difficult.  Furthermore,  because  of  variations  of  index  of 
refraction  in  the  atmosphere  and  other  outside  disturbances 
that  cause,  relative  motion  between  receiver  and  object,  active 
alignment  is  necessary.  Thus,  optical  trackers  are  needed  to 
keep  the  object  within  the  receiver's  field-of-view. 

Background 

Several  methods  are  currently  being  developed  in  the  Air 
Force  to  track  the  target  in  presence  of  several  disturbances. 

The  disturbances  include  target  motion,  mirror  vibration,  beam 
jitter,  stochastic  effects  of  the  atmosphere,  and  other  processes 
that  can  cause  relative  motion  between  the  beam  and  target. 

Method  of  Tracking.  One  scheme  to  accomplish  this  tracking 
task  uses  a direct  detection  optical  receiver  to  determine  the 
arrival  angle  of  the  incident  radiation  at  the  receiver.  The 
arrival  angle  is  converted  (through  optics)  into  a position  of 
a spot  of  light  on  the  active  surface  of  a photodetector  array. 

Or.ce  this  angle  arrival  is  determined,  it  can  be  used  to 
generate  error  signals  to  other  tracking  components  tn  correct 
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any  offset  alignment  between  the  receiver  and  target  that 
may  occur.  Because  of  disturbances  that  cause  relative  motion 
between  the  beam  and  target,  the  beam's  arrival  angle  varies 
at  the  optical  receiver.  Thus,  variations  of  arrival  angle  are 
therefore  converted  into  variations  of  the  spot's  position  on 
the  detector  surface.  The  spot  of  light  moves  dynamically  on 
the  surface,  releasing  photoelectrons  statistically  distributed 
about  the  spot's  position  of  maximum  intensity. 

A minimum -mean- square  estimator  which  will  be  known  as 
the  Snyder  Filter  estimates  the  position  of  maximum  intensity 
based  on  the  observations  of  the  photoelectrons  (Ref.  1).  These 
photoelectrons  are  measurement  data  to  the  filter  in  order  to 
obtain  an  estimate  of  the  spot's  current  position.  Here, 
the  Snyder  Filter  assumes  that  the  measurements  only  include  the 
target's  signal  intensity.  However,  other  disturbances  also 
cause  photoelectrons  to  be  released  from  the  surface.  These 
disturbances  which  degrade  the  Snyder  Filter's  tracking 
capability  are  caused  by  background  radiation  from  the  sky,  moon, 
sun,  stars,  and  other  sources,  and  dark  current  present  in  the 
absence  of  photoexcitation  on  the  photodetector.  Each  of  the 
noise  sources  are  discussed  in  the  next  chapter.  It  should  be 
noted,  however,  that  thermal  noise  is  ignored  for  photodetectors 
with  high  gain,  such  as  a photomultiplier  or  avalanche  photodiode. 

Quantum  Effects . At  higher  frequences,  the  quantum  nature 
of  light  is  an  important  consideration  in  many  optical  measure- 
ment systems  which  use  photon  or  quantum  detectors  to  detect 
these  optical  fields.  In  this  case,  light  can  be  described  as 
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a collection  of  discrete  packets  of  energy  called  photons. 

With  a given  probability,  each  photon  converts  into  an  electron. 
These  photoelectrons,  distributed  in  time  and  space,  are 
released  from  the  photodetector  surface.  Because  of  these  quantum 
effects  in  the  detection  process,  especially  when  the  field 
intensity  is  low  at  the  receiver's  input,  there  is  an  uncertainty 
associated  with  estimating  the  spot's  position  of  maximum 
intensity  based  on  the  sensed  photoelectrons.  Theoretically, 
both  the  quantum  nature  of  light  and  the  nature  of  the  detection 
process  fundamentally  limit  the  tracking  capability  of  optical 
trackers  (in  other  words,  the  ability  to  detect  the  signal 
photons  is  the  fundamental  limit). 

Because  of  the  quantum  effects  associated  in  the  detection 
process,  the  statistics  of  the  receiver's  output  process  must 
be  described  accurately;  the  Snyder  Filter  is  based  on  these 
statistics.  Experimental  evidence  indicates  that  the  statistics 
of  the  output  can  be  modelled  as  a random  point  process  in  time 
and  space  conditioned  on  knowing  the  optical  field  incident 
on  the  receiver  (Ref.  2).  A random  point  process  by  definition 
is  "a  mathematical  model  to  describe  a physical  phenomenon 

characterized  by  highly  localized  events  distributed  randomly 

u 

in  a continuum  (Ref.  3:2).  The  discrete  events  are  the  photo- 
electrons or  photoconversions  distributed  randomly  on  the  photo- 
detector surface.  Associated  with  each  event  are  the  temporal 
and  spatial  location  of  the  photoelectron.  Here,  the  event 
location  is  used  as  distributed  measurement  or  data. 
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in  order  to  estimate  the  current  position  of  maximum  intensity. 

To  summarize,  the  point  process  model  describes  the  fundamental 
uncertainty  of  photon -to-electron  conversion  in  the  measure- 
ment or  detection  process. 

Measurement  Models  (Gaussian  versus  Point  Process 
Descriptions).  To  contrast  the  point  process  measurement 
description  with  other  noise  models,  several  Gaussian 
descriptions  modelling  different  types  of  noise  are  briefly 
presented.  A stationary  white  noise  is  one  model  which  describes 
noise  as  having  a power  spectral  density  of  infinite  bandwidth, 
having  energy  at  all  frequencies.  Also,  infinite  bandwidth 
implies  that  the  noise  samples  are  uncorrelated  in  time. 

Although  no  noise  has  infinite  bandwidth,  the  white  noise  model 
does  provide  a reasonable  description  for  wideband  noise  having 
a wider  bandpass  than  the  physical  system  of  interest.  In  other 
words,  the  wideband  noise  appears  white  with  respect  to  the 
system.  In  many  instances,  however,  noises  are  not  well  modelled 
as  uncorrelated  in  time  (e.g.,  bandlimited  noise),  and  other 
models  must  be  exploited.  For  such  causes,  a Gaussian  white 
noise,  driving  the  input  of  a small  linear  system,  can 
generate  an  output  that  duplicates  or  closely  approximates  the 
second  order  characteristics  of  a given  noise  source.  The 

I linear  system,  or  "shaping  filter",  duplicates  the  second  order 

statistics  of  the  actual  noise  by  generating  a Gaussian  process 
with  the  same  second  order  statistics.  For  example,  the 
statistics  of  a zero-mean  exponentially  time-correlated  (First 
Order  Markov)  noise  can  be  duplicated  by  taking  the  output 
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statistics  of  a first  order  lag  driven  by  white-Gaussian  noise 
of  zero  mean.  This  model  approximates  a variety  of  experimentally 
observed  random  behavior  that  have  wide  or  narrow  bandwidths 
with  a flat  power  spectral  density.  Another  noise  model  is  the 
Second  Order  Markov  process  model  which  describes  phenomena 
having  vibration,  bending  or  any  other  periodic  random  charac- 
teristics. This  process  can  be  generated  as  the  output  of  a 
second  order  system  driven  by  white  noise. 

The  different  models  that  were  just  described  typically 
model  noise  as  due  to  a number  of  sources  that  are  always 
present  (or  continuous)  in  a dynamics  system  or  measurement 
device.  By  the  Central  Limit  Theorem  of  Statistics,  when  a 
number  of  random  variables  are  added  together,  the  sum  of  these 
random  variables  has  a probability  distribution  which  is  nearly 
Gaussian,  regardless  of  the  shape  of  the  individual  densities. 

The  Gaussian  models  are  exploited  in  many  cases  where  macro- 
scoptic  phenomenon  are  observed.  On  the  other  hand,  the 
quantum  effects  are  observed  as  a microscopic  interaction 
between  the  photons  in  the  incident  field  and  the  atoms  on  the 
photodetector  surface.  The  result  of  this  discrete  interaction 
is  the  probability  of  releasing  a photoelectron  at  a particular 
time  and  location  from  the  photodetector  array.  Thus,  various 
disturbances  affecting  a system  or  measurement  device  cannot 
always  be  described  as  having  additive  corrupted  noise  with 
Gaussian  statistics.  To  incorporate  the  quantum  effects  or 
uncertainty,  the  random  point  process  description  is  used  to 
model  the  measurement  process  from  the  optical  receiver  under 
limited  operations. 
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Assumptions  of  the  Snyder  Filter . This  report  assumes 
that  the  point  object  is  within  the  f ield-of- view  of  the 
optical  receiver  and  is  imaged  as  a single  spot  on  the  photo- 
detector  array.  In  other  words,  the  tracker  has  acquired  the 
target  and,  hence,  the  tracker  is  said  to  be  in  "fine-tracking" 
mode.  As  mentioned  previously,  the  solution  of  estimating  the 
position  of  maximum  intensity  based  on  released  photoelectrons 
is  the  Snyder  Filter.  For  this  research,  the  Snyder  Filter  is 
implemented  in  open-loop  configuration  (i.e.,  no  control  or 
feedback  input).  The  Snyder  Filter  models  the  spot  as  having 
a Gaussian-shaped  intensity  profile  and  the  spot  drifts 
dynamically  as  a First  Order  Gauss-Markov  process.  These 
characteristics  and  the  Snyder  Filter  are  discussed  in  detail 
in  the  next  chapter. 

In  addition  to  the  above  assumptions,  this  estimator 
assumes  its  measurements  are  from  a photodetector  array  capable 
of  measuring  the  "exact"  temporal  and  spatial  component  of  the 
photoelectron.  This  idealization  places  a limit  in  tracking 
capabilities  which  can  be  accomplished  with  photodetectors 
arrays  of  finite  resolution,  such  as  a quadrant  photomultiplier. 
However,  the  Snyder  Filter  assumes  that  its  measurement  data 
contain  no  extraneous  (or  noise)  measurements.  That  is,  photo- 
electrons are  not  released  due  to  background  radiation  or  dark 
current.  This  report  investigates  the  tracking  performance  of 
the  Snyder  Filter  based  on  observations  of  the  sensed  photo- 
electrons as  input  data. 
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Problem 


The  problem  is  to  establish  the  performance  capabilities 
of  the  Snyder  Filter.  In  this  case,  the  issue  is  whether  this 
minimum  mean-square- error  estimator  commits  errors  small  enough 
to  yield  good  tracking.  Thus,  the  problem  is  to  describe  the 
error  performance  of  the  Snyder  Filter  in  order  to  make  this 
evaluation.  Once  the  performance  description  is  established 
(under  the  assumptions  of  no  background  or  dark  current),  an 
additional  objective  is  to  determine  the  sensitivity  of  the 
Snyder  Filter  when  noise  data  is  presented  to  the  filter  as 
input.  To  improve  filter  performance  in  the  presence  of  the 
above  disturbances,  ad  hoc  methods  which  alter  the  filter 
structure  will  be  investigated. 

Overview  and  Approach 

Chapter  II  describes  the  models  employed  in  the  tracking 
scheme  and  establishes  conditions  in  tracking  performance  under 
which  the  Snyder  Filter  is  evaluated.  The  chapter  begins  by 
describing  the  model  of  the  direct  detection  optical  receiver 
and  its  detection  process.  In  this  section,  the  semiclassical 
approach  is  used  to  describe  the  photodetection  process.  The 
results  of  this  analysis  are  used  to  derive  the  statistics  of 
the  output  process  from  the  optical  receiver  (Appendix  A). 
Following  this  discussion  is  a description  of  the  detector 
noise  which  affects  the  output  of  the  receiver. 

Following  the  model  of  the  optical  receiver  is  a section 
on  the  Snyder  Filter.  A characterization  of  the  spot's 
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intensity  profile  and  its  dynamics  are  given.  Following  the 
subsection  is  a portrayal  of  position-sensitive  measurements 
available  to  the  Snyder  Filter.  As  noted  earlier,  these 
measurements  consist  of  point  process  photoconversions  in  time 
and  space.  Concluding  the  section  are  the  filter  equations 
of  the  Snyder  Filter. 

Chapter  II  finally  concludes  with  a description  of 
performance  parameters  to  describe  the  mean-square-error 
(MSE)  of  the  Snyder  Filter.  In  this  case,  performance 
descriptions  of  estimators  and  controllers  that  employ  point 
process  observations  in  time  and  space  are  available  in  recent 
literature  (Ref.  4).  However,  this  report  addresses  only  the 
estimation  performance  of  the  Snyder  Filter.  Also,  the 
Snyder  Filter  is  a minimum  mean-square-error  (MMSE)  estimator, 
under  ideal  conditions  (i.e.,  input  data  from  background 
radiation  or  dark  current  are  not  included).  The  MSE,  also 
denoted  as  EtP(t)l,  is  used  to  determine  whether  the  errors 
commited  by  the  Snyder  Filter  are  sufficiently  small  for 
tracking  purposes.  However,  the  MSE  cannot  be  evaluated  in 
closed  form;  thus,  upper  and  lower  bounds  are  derived  rather 
than  solving  for  MSE.  In  Appendix  B,  the  models  and  conditions 
of  the  Snyder  Filter  are  applied  in  the  upper  and  lower  bounds. 

In  the  appendix,  two  physically  motivated  parameters  are 
defined  from  these  bounds  to  establish  conditions  under  which 
the  Snyder  Filter  yeilds  good  tracking:  the  average  number 
of  photoelectrons  in  a coherence  time  of  the  spot's 
dynamics  and  the  mean  square  value  of  the  spot's  jitter 
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normalized  to  the  square  of  the  beamwidth.  These  parameters 
establish  regions  of  operations  where  the  upper  and  lower 
bounds  converge  to  the  actual  MSE  (i.e.,  bounds  are  identical). 

The  parameters  are  discussed  in  detail  at  the  end  of 
Chapter  II.  Note  that  the  performance  measures  themselves 
assume  perfect  data;  that  Is,  these  measures  of  performance 
are  not  valid  for  cases  which  include  data  from  dark  current 
and  background  radiation. 

To  demonstrate  the  tracking  performance  provided  by  the 
above  parameters  and  to  determine  the  filter's  sensitivity 
to  noise  measurements,  a Monte  Carlo  simulation  and  performance 
analysis  is  performed  on  a digital  computer.  Chapter  III  develops 
algorithms  to  simulate  the  models  in  the  computer  required  to 
perform  this  study.  The  chapter  begins  with  a section  on  simu- 
lating the  spot's  dynamics  on  the  detector  surface  and  the 
required  noise  generators  are  discussed  next.  Following  these 
sections  are  time-space  point  process  generators.  Two  methods 
are  used  to  generate  these  data  measurements:  one,  the  Binomial 
Method,  uses  the  incremental  definition  of  the  point-process, 
and  the  other,  the  Conditional  Method,  involves  Bayes'  rule 
and  the  statistics  of  the  output  of  the  optical  receiver.  The 
Binomial  Method  involves  a threshold  test  to  declare  an  event 
(.photoconversion)  at  each  infinitesimal  area  on  the  detector 
surface.  The  Conditional  Method  calls  a uniform  random  variable 
distributed  between  zero  and  one  from  the  computer  system  and 
transforms  it  into  another  random  variable  with  a desired 
probability  distribution.  This  transformed  variable  is  a 


10 


realization  of  either  the  time  or  space  coordinates  associated 
with  the  photoconversion  process.  After  describing  these  two 
methods,  the  theoretical  and  practical  aspects  of  implementing 
each  method  in  the  computer  are  discussed.  Concluding 
Chapter  III  is  the  performance  analysis  which  characterizes  the 
error  process  generated  by  the  Snyder  Filter.  Given  in  this 
section  are  the  sample  statistics,  plots,  and  test  cases  required 
to  evaluate  filter  performance.  Assuming  perfect  measurements 
(i.e.,  no  background  or  dark  current  events),  the  actual  per- 
formance of  the  Snyder  Filter  is  compared  with  the  upper  and 
lower  bounds  on  MSE.  This  analysis  should  indicate  where  the 
actual  performance  lies  with  respect  to  the  bounds.  Then 
noise  measurements  are  included  in  the  analysis  to  determine 
the  filter  sensitivity  to  these  measurements.  Finally,  ad  hoc 
methods  of  residual  monitoring  and  filter  tuning  are  employed 
to  enhance  tracking  performance  in  the  presence  of  the  noise 
measurements. 

Chapter  IV  contains  results  from  the  simulation.  For 
perfect  measurements,  the  simulation  results  indicate  the  useful- 
ness of  the  bounds  in  describing  the  actual  performance  of  the 
Snyder  Filter.  Also,  these  results  show  for  low  signal  events 
per  coherence  time  that  the  upper  bound  describes  the  actual 
performance.  When  measurements  from  background  radiation  or 
dark  current  are  included  in  the  simulation,  the  estimator 
was  found  to  be  quite  sensitive  to  these  noise  measurements, 
resulting  in  very  poor  tracking.  Finally^  it  was  found  that 


tracking  performance  can  be  improved  substantially  by  an 
ad  hoc  procedure  called  residual  monitoring. 

The  last  chapter  is  a summary  of  conclusions  and 
recommendations  for  further  research. 


II.  The  Snyder  Filter  and  Its  Open  Loop 
Estimation  Performance 

The  Snyder  Filter  is  a recursive  data  processing 
algorithm.  This  algorithm  employs  point  process  observations 
in  time  and  space  to  estimate  optimally  the  position  of  maximum 
intensity  on  a detector  array  of  infinitesimally  fine  partition. 
However,  its  estimate  is  only  as  good  as  the  model  describing 
the  output  statistics  of  the  photodetector  array  (the  key 
component  in  the  optical  receiver).  Therefore,  the  model  of 
the  direct  detection  optical  receiver  and  its  detection 
process  deserves  description.  The  modelling  will  involve  a 
statistical  description  of  the  receiver's  output.  Following 
the  receiver  model,  the  Snyder  Filter  is  discussed.  Concluding 
the  chapter  is  the  tracking  performance  by  which  this  estimator 
is  to  be  judged  using  upper  and  lower  bounds  on  MSE.  For  this 
case,  two  parameters  are  defined  to  describe  the  actual  per- 
formance of  the  Snyder  Filter.  The  concepts  are  presented 
in  one  spatial  dimension. 

Optical  Receiver  Model 

There  are  two  classes  of  optical  sceivers,  direct 
detection  receivers  and  heterodyning  i ”.  elvers.  The  difference 
between  these  two  receivers  is  that  the  former  operates  to 
detect  the  instantaneous  power  in  the  collected  field  at 
receiver's  input  and  the  latter  operates  to  detect  a combined 
field  formed  by  mixing  a locally  generated  field  with  the 
received  field.  This  study  is  concerned  with  direct 
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detection  receivers. 


General  Description.  The  direct  detection  receiver  is 
also  known  as  a noncoherent  detector  which  is  shown  in  Fig. 

2.  The  receiver  consists  of  a lens  system  and  a photo- 
detector. The  lens  system,  also  referred  to  as  the  front 
end  of  the  receiver,  performs  two  functions.  One  is  to 
focus  the  optical  field  onto  the  photodetector  surface  at 
a focal  point.  Another  function  is  to  provide  some  degree  of 
reducing  the  intensity  of  the  background  noise  (i.e.,  radia- 
tion from  the  sun  or  sky)  prior  to  photodetection.  The  photo- 
detector then  converts  the  incident  field  into  an  electrical 
signal. 

Photodetectors . Photodetectors  consist  of  two  types: 
photon  and  thermal  detectors.  Because  thermal  detectors  have 
a slower  response  for  the  frequencies  of  interest,  they  are 
not  considered  further  (Ref.  5:87).  On  the  other  hand,  photon 
or  quantom  detectors  respond  quickly  at  the  optical  frequencies. 
An  excellent  discussion  on  quantum  detectors,  their  description, 
performance,  and  characteristics,  is  given  in  Refs.  5:87  and 
6:298  . However,  of  all  the  available  types  of  detectors  using 
photosensitive  materials,  all  behave  according  to  quantum- 
mechanical  principles  (Ref.  2:6).  Figure  3 portrays  the  photo- 
detector model.  These  detectors  use  photosensitive  materials 
to  produce  electrical  signals  in  response  to  changes  in  the 
input  field  intensity.  The  signals  are  electrons,  released 
from  the  photodetector  surface,  and  they  are  collected  at 
the  anode  by  an  electric  field.  This  collection  and  movement 
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of  charge  at  the  anode  creates  an  output  current.  Because 
the  conversion  of  optical  field  to  electron  flow  and  the 
arrival  times  of  electrons  at  the  anode  are  probabilistic,  the 
observed  output  current  is  a random  process  in  time  and  space 
(Ref.  2:39).  This  fact  is  true  even  if  the  input  field  is 
deterministic.  In  short,  the  detection  operation  inherently 
induces  a randomness  in  the  detector  output. 

To  model  the  relationship  between  the  optical  field  and 
the  number  of  photoelectrons  released  involves  a treatment 
in  quantum  mechanics  (Ref.  3).  This  treatment  has  two 
approaches  to  describe  the  above  phenomenon.  One  is  to 
describe  the  field  as  a collection  of  discrete  packets  of 
energy  called  photons.  Here,  each  photon  produces  an  electron 


15 


Figure  3.  Photodetector  Model 

(copied  from  Ref.  2:38) 

at  the  detector  surface  with  a given  probability.  This 
probability  of  conversion  is  termed  the  quantum  efficiency, 
n . The  second  method,  referred  to  as  a semiclassical  analysis 
(a  direct  result  of  quantum-mechanical  considerations),  treats 
the  focussed  field  as  a wave, and  it  incorporates  a statistical 
relation  to  describe  the  interaction  between  the  optical  wave 
and  the  detector  surface.  The  second  method  is  preferred  for 
its  simplicity  and  the  added  insight  it  provides.  An  outline 
of  the  second  approach  is  found  in  Ref.  2 and  a detailed 
description  is  given  in  Ref.  7 . Since  a detailed  description 
of  the  semiclassical  analysis  is  quite  lengthy,  only  results 
of  this  semiclassical  approach  are  given.  These  results  are' 
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as  follows:  in  a short  infinitesimal  time  interval.  At, 
the  release  of  an  electron  from  an  atom  at  a differential 
lenght  AX.  on  the  detector  is  proprotional  to  the  field 
intensity  over  the  observation  area,  AX At.  That  is. 


‘one  electron  released  * 
in  a length  AX  during 
the  time  At 


- y!(t,r)AtAX 


AXAt+O 

where  P[-]  denotes  probability  of,  y is  a proportionality 
constant,  I(t,r)  denotes  the  field  intensity  at  time  t and 
point  r on  the  detector  surface,  and  AX  is  located  at  r. 
Other  results  from  this  approach  are  as  follows: 


no  electrons  released 
from  AX  during  At 


1-y I (t , r ) At AX  (2) 


AXAt-K) 


[more  than  one  electron" 
released  from  AX 
during  At 


AXAt->0 


Along  with  these  results  is  an  assumption  that  the  release  of 
electrons  from  disjoint  observation  areas  At AX  are  treated 
as  independent  events.  This  assumption  and  Eqs  (1)  through 
(3)  mathematically  model  the  photodetecting  surface. 

Detector  Noise . With  no  noise  sources,  the  output 
current  is  given  by 


i(t,r)  = q£  6 (t-t,  ) 6 (r-r,  ) 
k=l  K K 
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t h 

where  q is  the  electronic  charge,  t^  is  the  ku  event 
arrival  time,  r^  is  the  corresponding  event  location,  r 
denotes  a position  on  the  detector  surface,  n is  the  number 
of  events  in  an  observation  area,  and  6 ( t ) and  <S(r)  are 
delta  functions.  Eq.  (4)  assumes  that  the  detector  has  infinite 
bandwidth  (e.g.,  a photomultiplier  or  avalanche  photodiode 
fits  closely  to  this  description).  The  relation  also  indicates 
that  the  current  is  a random  process  since  n,  t.  , and  r.  are 
random  variables. 

However,  the  imperfection  of  the  optical  receiver  and 
background  radiation  also  affect  the  output  current.  In 
general,  the  output  current  from  the  photodetector  array 
at  (t,r)  can  be  represented  as 

i(t,r)  = ig(t,r)  + ib(t,r)  + id(t,r)  + ith(t,r)  (5) 

where  ig(t,r)  is  due  to  signal  current,  ib(t,r)  is  current 
due  to  background  noise,  id(t,r)  is  due  to  dark  current,  and 
ith(t,r)  is  current  due  to  thermal  noise.  The  dark  current 
is  a result  of  current  flowing  in  absence  of  photoexcitation. 
Also,  the  dark  current  is  well  modelled  as  a point  process 
(Ref.  5:148)  with  a rate  function  that  is  constant  (i.e., 

Xd(t,r)  = is  termed  homogeneous (Ref . 3:54)).  Further, 

the  background  noise  can  be  associated  with  a point  process 
in  nighttime  operation,  and  its  rate  function  is  homogeneous 
provided  that  three  conditions  are  satisfied:  (1)  background 
electric  field  has  zero  mean,  (2)  background  intensity  is 
constant,  and  (3)  the  product  between  the  ensemble  average 
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of  the  number  of  electrons,  n^,  and  the  bandwidth  of  the  back- 
ground radiation,  — , is  much  less  than  1.  The  last  condition, 

Tb 

n^Tj3<<l,  is  known  as  low  photon  coherence  which  is  a valid 
approximation  for  nighttime  operation.  For  convenience,  the 
dark  current  in  this  study  will  be  the  combined  effects  due 
to  the  detector  dark  current  and  due  to  background  radiation. 
Finally,  the  thermal  or  Johnson  noise  current  i^^(t,r)  is 
caused  by  thermal  fluctuations  of  electrons  in  a resistor. 

It  is  modelled  as  a white  Gaussian  random  process  (Ref.  5:1^5). 

To  describe  the  output  current  in  Eq.  (5)  statistically 
is  very  difficult.  It  is  difficult  because  the  output  current 
is  a mixed  process  (i.e.,  a sum  of  a discrete  process  (Point- 
Process  model)  and  a continuous  process  (Gaussian  model)). 
Hence,  limitirg  cases  must  be  considered.  This  study  will 
ignore  the  current  due  to  thermal  noise  since  one  of  the  objec- 
tives of  this  study  is  to  observe  the  quantum  effects  of  dark 
current  in  filter  performance.  Also,  if  the  detector  is  a 
photomultiplier  or  avalanche  photodiode  with  high  gain,  one 
can  assume  for  low  input  intensity  that  the  optical  receiver 
is  operating  in  the  point  process  regime,  and  thus,  thermal 
noise  is  ignored. 

Photodetector  Statistics . It  is  assumed  In  the  detector 
model  that  the  electrons  arriving  at  the  anode  of  each  detector 
of  an  array  are  independent.  Then  the  output  statistics  of 
each  detector  can  be  treated  independently  of  one  another. 

One  needs  the  statistics  describing  the  number  of  events,  n, 
the  ordered  event  times,  t^'s,  and^ordered  event  locations. 
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r^'s,  to  characterize  the  detector  output  completely.  Thus, 
the  joint  density  of  these  variables  is  needed.  With  the 
results  given  in  Eqs.  (1)  through  (3),  the  above  joint 
density  of  the  photodetector  can  be  derived.  Before  stating 
the  results,  it  is  convenient  to  define  a complex  envelope 
which  characterizes  (in  part)  the  focussed  or  received  field 
(i.e.,  a classical  description).  Let  u(t,r)  denote  the  scalar 
optical  field.  Then  the  complex  envelope  U(t,r).of  the  scalar 
field  is  defined  implicitly  as 

u(t,r)  = Re(U(t,r)eJ'2irf«t)  (6) 

where  R(*)  denotes  the  "real  part  of"  the  quantity  in  paren- 
thesis, and  fg  denotes  the  carrier  frequency.  The  intensity 
of  the  field  is  represented  as 

Kt.r)  = \ \ |U(t,r)  |2  (7) 

m 

where  | * | denotes  the  modulus  of  the  complex  field,  and  Z 
denotes  the  impedance  of  the  medium  (Ref.  2:10). 

With  the  definition  of  a complex  envelope  and  using 
Eqs.  (6)  and  (7),  the  output  of  the  photodetector,  conditioned 
on  knowing  the  optical  field,  can  be  modelled  as  a point  process 
(Ref.  2)  with  rate  function  (in  one  dimension) 

X(t,r)  = j z\  U (t,r)|2di  + (8) 

due  to  signal  due  to  dark  current 

(Xd(t,r)  = Xd) 
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where  n is  quantum  efficiency,  h is  Planck's  constant,  f is 
the  optical  signal  frequency,  U (t,r)  is  the  complex  envelope 
of  the  scalar  field  incident  on  the  detector  surface  at  (t,r), 
l is  the  length  of  detector  array,  and  di  is  the  differential 
detector  length.  The  integral  in  Eq.  (8)  represents  the 
total  optical  power  deposited  by  the  beam  on  the  detector 
surface  and  the  energy  of  a photon  is  hf  . Also,  recall  from 
the  previous  section  that  the  homogeneous  dark  rate  function 
is  due  to  the  combined  effects  of  background  radiation  (of 
zero  mean  electric  field)  satisfying  low  photon  coherence  and 
dark  current  resulting  from  imperfection  of  the  photodetector. 
This  dark  rate  function,  denoted  as  A^,  which  is  assumed  to 
be  independent  in  time  and  space  can  be  expressed  as 


due  to  background  due  to 

radiation  photodetector 


Note  that  in  Eq.  (8)  and  (9)  the  optical  receiver  assumes  to 
be  operating  at  the  point-process  regime.  Also,  the  first 
term  at  the  right-hand  side  of  Eq.  (9)  is  a constant  rate 
function,  since  this  study  assumes  a constant  background 
intensity.  In  other  words,  the  dark  current,  A^,  is 
independent  of  time  and  space.  Finally  with  Eqs . (1)  through 
(3),  the  output  statistics  of  the  optical  receiver  conditioned 
on  knowing  the  input  field  is  found  in  Appendix  A as 
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f[{tk},{rk},  n]  = exp ( -/£/q  ^(a.B)dctde)  n X(tk,rk)  (10) 

k— 1 


where  {tk)  denotes  set  of  ordered  event  times,  {rk>  denotes  the 
corresponding  set  of  event  locations,  1 is  detector  length, 
n is  the  number  of  events  in  a time  interval  [0,T]  on  £.  Again, 
this  joint  density  completely  characterizes  the  statistical 
output  of  the  optical  receiver. 

The  Snyder  Filter 

The  previous  sections  described  the  statistical  model 
of  the  optical  receiver  which  is  used  in  the  optical  tracking 
problem.  The  lens  system  of  a receiver  focusses  the  received 
field  on  a detector  surface  at  the  focal  plane.  The  array  of 
detectors  responds  to  the  received  field  which  arrives  at 
various  angles  at  the  receiver.  Because  of  random  disturbances, 
the  spot  wanders  randomly  on  the  detector  plane.  Electrons 
are  released  at  a rate  proportional  to  the  light  intensity, 
given  by  Eq.  (8).  For  the  model  upon  which  the  Snyder  Filter 
is  based,  xd  is  assumed  to  equal  zero.  The  Snyder  Filter 
estimates  the  position  of  maximum  intensity  on  the  detector 
plane  in  terms  of  these  time-space  photoelectrons.  An  outline 
of  its  models  and  results  will  now  follow  specializing  for 
one  spatial  dimension. 

Characterization  of  Spot . The  spot  of  light  on  the 
detector  surface  is  a diffraction  pattern  focussed  by  the 
receiver's  lens.  It  is  assumed  that  the  spot  of  light  has 
a Gaussian-intensity  profile  given  by 
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I(r,x(t),t)  = I o ( t ) exp[-g(r2^^  | )2]  (li; 

where  Ig(t)  is  the  maximum  intensity,  o(t)  characterizes  the 
beam  spread,  and  x(t)  is  the  true  location  of  the  spot 
(i.e.,  the  location  at  the  centroid  of  the  beam).  For  this 
study,  I q C t ) and  o(t)  are  assumed  to  be  static  (i.e.,  Ig(t)  = 
IQ  and  o(t)  = o).  Figure  4 portrays  this  intensity  profile. 
From  the  figure,  the  beam  width  (or  beam  radius)  in  this 


study  is  defined  as  a point  from  the  maximum  intensity,  IQ, 

-1/2 

to  a point  in  which  the  intensity  is  reduced  to  Ine 

* 0 , or 


about  60.7%  of  the  maximum  intensity. 

Characterization  of  Spot ' s Motion.  The  motion  of  the 
spot's  centroid  is  assumed  to  be  modelled  by  the  following 
linear  stochastic  differential  equation 


dx(t)  = F(t  )x(t ) + G(.t)dg(t);  x(tQ)  = xQ 


where  §(t),  t>0  is  a standard  Wiener  process,  F(t)  is  the 
state  coefficient,  and  G(t)  is  the  noise  input  coefficient. 

The  spot's  dynamics  are  governed  by  a number  of  disturbances. 
The  atmosphere  is  one  such  disturbance;  characterized  as  an 
inhomogeneous  random  optical  channel  (Ref.  9:16-26),  the 
atmosphere  changes  its  index  of  refraction  due  to  temperature 
variation.  Because  of  the  stochastic  nature  of  the  atmosphere, 
and  frequency  at  which  the  field  propagates,  the  optical  wave 
interacts  with  the  foreign  particles  in  the  atmosphere.  One 
consequence  of  this  phenomenon  is  that  the  received  field 
causes  variations  in  the  beam's  arrival  angle,  causing  the 
spot  to  move  randomly  on  the  detector  plane.  Platform 
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vibration  of  the  optical  receiver  is  another  disturbance  which 
causes  the  spot  to  drift.  The  important  point  is  that  any 
relative  motion  between  object  and  receiver  due  to  other 
effects  will  result  in  the  spot  to  drift  dynamically  on  the 
detector  plane.  However,  it  is  not  the  purpose  of  this  study 
to  model  accurately  the  motion  of  the  spot,  but  to  provide  a 
reasonable  and  simple  model  for  digital  simulation.  A model 
was  developed  to  describe  the  dynamics  of  the  target  relative 
to  the  receiver  (Ref.  10:42).  This  model  is  a stationary 
First  Order  Gauss-Markov  process,  the  output  of  a first  order 
lag  driven  by  white  Gaussian  noise  of  zero  mean.  This  process 
also  describes  a variety  of  phenomena  having  bandlimited 
characteristics.  Thus,  this  model  will  be  used  to  describe  the 


dynamics  of  the  spot.  Another  name  for  this  process  is  an 

exponentially  time-correlated  Gaussian  process.  Figure  5 

portrays  the  properties  of  a First  Order  Markov  Process.  Two 

parameters  of  interest  to  describe  the  output  of  the  first 

order  lag  are  the  coherence  time  x , and  its  mean  squared 
2 

value  ox  . The  coherence  time  (or  correlation  time),  describes 
how  fast  the  resulting  output  process,  x(t),  will  vary  in  time. 
The  reciprocal  of  tc  is  viewed  as  the  bandwidth  of  the  distur- 
bance process  (i.e.,  turbulance  plus  other  effects).  Heuris- 
tically,  the  spot's  rms  value  of  random  displacements  (or 
spot's  rms  jitter),  is  a measure  of  how  much  the  spot  will 
fluctuate  about  its  mean  position  on  the  detector  array. 

Finally,  the  strength  of  a white  Gaussian  noise  driving  the 
first  order  lag  is  given  as 

E[w(t)w(t+x)]  = Q5(x)  (13) 

where  Q = — a 2 (Ref.  9:4.77).  The  strength,  Q,  describes 
Tc  x 

how  fast  the  Wiener  Process  or  Brownian  Motion,  6(t),  diverges 
in  mean  square  value  from  its  initial  condition  of  zero  (Ref. 
9:4.29) . 

Position-Sensitive  Measurements . As  mentioned  earlier  in 
the  report,  these  measurements  consist  of  event  locations 
which  provide  information  about  the  position  of  maximum 
intensity.  Assuming  no  dark  current,  the  conversion  rate 
at  the  detector  surface  is  modelled  as 
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Figure  5.  Properties  of  a First  Order  Markov  Process 
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A (r, x(t ) ) 
g ~ 


_/nIo/27°\r  1 exp{^(r-x(t))2}]  (14) 

\ hfQ  /L/27o  20  J 

If  the  detector  length  is  sufficiently  large,  the  quantity 
in  brackets,  [ - J , in  Eq.  (14)  integrates  over  the  entire 
detector  length  to  a value  very  close  to  one.  Note  that  Eq. 

( 1 4 ) is  derived  by  substituting  Eq.  (11)  as  the  signal  intensity 
in  the  first  term  of  Eq.  (8).  The  term  ri Ig /2tF  a,  represents 
the  average  optical  energy  deposited  on  the  detector  by  the 
spot  to  release  photoelectrons.  The  term  hfg  represents  the 
energy  of  a single  photon,  characterizing  the  particle 
nature  of  light  of  frequency  f q . Then  integrating  Eq.  (14) 
both  over  time  (with  limits  from  0 to  T)  and  space  (over 
entire  detector  length),  yields  (77^-  /2tT  a )T.  This  term 

nig 

represents  the  ensemble  average  of  the  number  of  photons 
detected  at  the  output  over  the  time-space  observation  area. 
Recall  that  the  data  measurements  are  obtained  from  a photo- 
detector array,  capable  of  measuring  the  exact  temporal  and 
spatial  component  of  a photoelectron.  With  this  idealization 
the  estimator’s  performance,  derived  in  Appendix  B,  places  a 
limit  in  performance  which  can  be  achieved  with  arrays  of 
finite  resolution.  In  other  words,  a decrease  in  spatial 
resolution  should  result  in  poorer  performance.  It  is  also 
assumed  that  the  detector  array  is  large  such  that  edge  effects 
of  the  spot's  intensity  profile,  described  as  Gaussian-Shaped 
can  be  ignored.  A good  "rule  of  thumb"  is  to  have  the  detector 
large  such  that  the  distance  between  the  center  of  the  spot 
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and  the  edge  of  the  detector  array  is  at  least  six  beam 
widths.  This  is  "good"  because  with  respect  to  a Gaussian- 
shaped  beam,  the  above  "rule  of  thumb"  closely  approximates 
a detector  of  infinite  length.  Also,  by  having  the  detector 
length  large#  edge  effects  may  be  ignored  as  the  spot  moves 
dynamically  even  with  large  wanderings  on  the  detector  surface. 

Snyder ' s Filter  Equations . This  subsection  presents  the 
dynamical  filter  equations  that  incorporate  point-process 
observations  to  estimate  the  position  of  maximum  intensity. 

These  equations,  specialized  for  one-dimension  state  as  follows: 

dx(t)  = F(t)x(t)dt  + /^K  (t ) [r-H( t ) x (t ) N(.dtxd£)  (15) 

dP(.t)  = F(t  )P  (t  )dt  + P(t)F(t)dt  + G(t)Q(t)G(t)dt 

- /AK(t)H(.t)N(dtxdA)  (16) 

K(t ) = P(t)H(t)[H(.t)P(t)H(.t)  + R(.t)]_1  C17) 

x(tQ)  = 0,  PCtQ)  = Pg  (initial  conditions  at  tQ)(l8) 

where  x(t)  is  the  conditional  estimate  of  the  spot's  centroid 
given  the  measurement  data  Ci.e.,  event  locations  that  are 
random  in  time  and  space),  P(t)  is  the  conditional  error 
variance,  H(t)  is  the  measurement  coefficient,  K(t)  is  the 
filter  gain,  r is  the  location  of  the  photoconversion  from  a 
differential  detector  length  d l in  the  differential  time, 
dt,  and  N(dtxd£)  = dN  is  the  differential  number  of  events 
in  the  differential  "space-time"  observation  area,  dtxdJl. 
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The  filter  equations,  Eq.  (15)  through  ( 1 8 ) are  recursive 

relations,  similar  in  form  to  the  discrete-time  measurement 

Kalman-Bucy  Filter.  They  differ  in  that  the  Snyder  Filter 

updates  its  estimate  and  uncertainty  at  the  event  times 

whereas  the  Kalman-Bucy  Filter  updates  at  predetermined  times. 

The  mean  and  the  variance  of  the  conditional  probability 

density  function  which  describes  the  error  in  the  estimate 

of  the  spot's  centroid  are  propagated  forward  from  the  time 
t h 4" 

period  of  the  k photoconversion  tk  , to  the  time  of  the 
(k+l)“st  photoconversion,  tk+1~.  Then  the  filter  updates  at 
time  tk+^,  by  processing  the  measurement  at  that  time,  to  achieve 

■*  4. 

the  estimate  x(tk+^  ),  where  denotes  before  measurement 
update  and  "+"  denotes  after  measurement  update 

It  is  convenient  to  expand  further  the  interpretation 
of  the  filter  equations.  During  interarrival  times  between 
photoconversions,  the  filter's  conditional  estimate  of  the 
centroid's  position  x(t),  and  Its  error  variance  P(t),  evolve 
in  time  by  the  following  equations 

^r-  = nt)xLt)  ; tk<t  < tk+1  (19) 

§f  = 2F(t)P(t)  + G(t)Q(t)G(t)  ; tfc<t  < tk+1  (20) 

The  above  relations  govern  how  x(t)  and  its  associated  conditional 
error  variance  P(t)  will  propagate  when  the  filter  does  not 
observe  a photoconversion.  But  when  the  filter  dees  observe  an 
event,  via  the  detector  array,  it  will  update  x(t)  and  P(t) 
after  the  occurrence  of  the  event  by  the  following  update 
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equations : 

P(tk+)  = [H(tk)R'1(tk)H(tk)  + P(tk“)J_1  (21) 

x(tk+)  = P(tk)[H(tk  )R_1(tk)rk  + P_1(tk")x(tk")]  (22) 

where  tk  is  the  time  after  the  occurrence  of  the  k photocon- 
version, tk~  is  the  time  before  this  occurrence , and  rk  is 
the  location  of  photoconversion  Cor  photoelectron).  By  using 
the  matrix  inversion  lemma  (Ref.  1),  Eqs  (21)  and  (22)  may  be 
alternatively  expressed  as 

x(tk+)  = x(tk")  + K(tk“) [rk-H(tk)x(tk-)]  (23) 


where 


P(tk+)  = P(tk")  - K(tk~)H(tk)P(tk")  (24) 

K(tk)  = P(tk“)HCtk)[H(tk)P(tk-)H(tk)  + R(tk)]  1 (25) 


Eqs.  (22)  and  (23)  indicate,  clearly,  the  recursive  nature  of 

A A 

how  x(t)  and  P(t)  evolve  in  time  given  the  history  of  data  ' 
measurements.  In  other  words,  the  "recursive"  nature  of  the 
Snyder  Filter  does  not  require  all  previous  data  measurements 
to  be  kept  in  storage  and  processed  every  time  a photoconversion 
is  observed.  This  saves  computer  resources  which  is  an  impor- 
tant consideration  for  filter  implementaton. 

To  summarize,  the  Snyder  Filter  is  a recursive  data 
processing  algorithm,  using  the  available  information  from 
the  photoconversions  to  improve  the  accuracy  of  estimating 
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the  spot's  current  position  of  maximum  intensity.  The  Snyder 
Filter  is  a minimum  mean  square  error  estimator;  it  attempts 
to  provide  the  best  estimate  possible  in  such  a manner  that 
the  error  is  minimized  statistically. 

The  next  section  specifies  parameters  where  sufficient 
conditions  are  presented  to  provide  actual  mean-square-error 
(MSE) . 

c I 

Upper  and  Lower  Bounds 

A question  should  be  asked  whether  or  not  the  filter 
errors  are  small  enough  for  tracking  purposes.  Thus,  the 
MSE  (i.e.,  an  expectation  of  the  conditional  error  variance 
given  the  data  measurements  over  all  possible  histories  of 
measurements  taken),  is  used  to  make  this  evaluation. 

Unfortunately,  finding  the  mean-square-error  cannot  be  evaluated 
in  closed  form.  Rather  than  attempt  to  solve  for  the  MSE, 
one  can  bound  it  above  and  below.  The  method  involves  easily 
precomputed  upper  and  lower  bounds  (P#4  lower  bound  on  MSE  and 
P*A  upper  bound  on  MSE). 

The  following  is  a quote  from  Ref.  (4)  in  how  these 
bounds  are  derived: 

"The  upper  bounds  are  derived  by  evaluating 
exactly  the  performance  of  a parameterized  family 
of  suboptimum  design;  one  of  these  is  identified 
as  having  smaller  performance  than  any  other,  thus 
providing  a minimal  upper  bound  within  this  family 
of  suboptimum  design.  The  lower  bounds  are  obtained 
directly  by  calculations  involving  inequalities." 

The  above  results  are  specialized  for  this  open-loop  problem 

in  one  dimension.  Here,  the  upper  and  lower  bounds  are 
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solved  In  steady  state  (i.e.,  = 0)  and  are 

derived  in  Appendix  B.  The  solution  to  this  steady-state 
estimation  problem  involves  a quadratic  equation  yielding  a 
result  given  as 


# 

~2 


* 

P 

~2 

a 


2X 


£-1  + (1  + ^4XY ) 1/2  J 


(26) 


= 2(Y+1)'  [(Y-X)  + ( (Y-l ) 2 + 4(X+1)Y}1/2J  (27) 


where 


X = 


Y = 


(h^  ^ o)(1/2tc) 


(28) 


(29) 


Several  comments  deserve  mentioning  with  regard  to  Eqs. 

(26)  through  (29).  The  upper  and  lower  bounds  are  normalized 
with  respect  to  the  square  of  the  beam  width.  This  normalization 
is  convenient  since  the  error  in  estimating  the  spot  car 
be  measured  in  terms  of  beam  width  (i.e.,  a performance  or 
criterion  can  be  established).  Thus,  for  good  tracking 
performance,  it  is  desirable  to  have  and*--^^  small. 

Figures  (6a-6i)  are  plots  of  the  upper  and  lower  bounds  on 
MSE  as  a function  of  X and  Y.  Note  the  absolute  magnitude 
and  relative  separation  of  the  bounds  as  X and  Y are  varied 
on  each  of  these  plots.  Consider  Figs.  6a-6d  : for  a fixed 
Y,  the  bounds  both  decrease  in  magnitude  and  in  relative 
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Figure  6d.  Upper  and  Lower  Bounds  on  MSE 
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seperation  as  X gets  larger.  Physically,  the  X-pararaeter  can 
be  described  as  the  number  of  photoconversions  in  a coherence 
time  of  the  spot's  dynamics.  As  the  number  of  data  observa- 
tions increases  (i.e.,  X-parameter  increases),  the  filter 
has  more  information  to  correct  and  update  its  estimate  of 
the  spot's  current  position  of  maximum  intensity.  Thus,  the 

filter  can  update  its  estimate  frequently  if  a laf^e  number 

l 

of  events  are  observed  in  a relatively  short  period  of  time. 

Suppose,  for  example,  the  spot  drifts  on  the  detector  with 

large  "wanderings".  Also,  assume  the  filter  has  just  processed 

one  event  and  does  not  observe  another  event  in  a relatively 

long  time  (i.e.,  a low  rate  function  or  a small  X).  Then  the 

information  provided  by  the  "processed"  event,  has  little 

value  at  some  later  time  if  the  spot  is  "jittering"  a; 

higher  frequency  with  large  "wanderings"  (i.e.,  spot 

dynamics  have  a small  coherence  time,  x , and  a large  jitter, 

2 

o ).  On  the  other  hand,  if  the  filter  frequently  observes  a 

A 

larger  number  of  events  within  a coherence  time,  the  Snyder 
Filter  has  more  information  concerning  the  spot's  location  on 
the  detector  surface. 

Consider  now  a fixed  X and  a variation  of  the  Y-parameter 
(Fig.  6e-6i).  As  Y increases,  the  bounds  increase  in  magnitude 
and  in  relative  separation.  Y is  the  ratio  between  the  spot's 
jitter  normalized  by  the  square  of  the  beam  width.  The  beam 
width  indicates  whether  the  filter  should  track  the  event 
locations  closely  in  order  to  locate  the  spot.  As  the  beam 
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width  increases,  photoelectrons  are  distributed  at  larger 
distances  from  the  spot's  centroid.  Hence,  increasing  the 
beam  width  results  in  having  less  weight  in  the  point-process 


measurements.  Thus,  the  optimal  gain,  K(t),  should  decrease 
as  seen  in  Eq.  (17). 

2 

Similarly,  the  noise  strength,  ox  , indicates  stronger 

noise  sources  driving  the  dynamics  of  the  spot  or  increasing 

2 

the  uncertainty  in  the  model.  Increasing  ox  increases  both 
the  rate  of  growth  of  the  conditional  error  covariance,  P(t), 
between  event  times.  This  results  in  increasing  values  of 
the  filter  gain,  implying  larger  weighting  of  the  point  process 
measurements.  Stated  another  way,  increasing  t'fte  jitter 
indicates  that  the  filter  must  put  less  "confidence"  in  the 
prediction  capability  of  the  filter's  own  internal  dynamics 
model. 

I 

Physically,  it  is  desirable  to  have  a narrow  beam  width 
such  that  electrons  can  be  distributed  near  the  position  of 
maximum  intensity.  At  the  same  time,  the  spot’s  jitter 
should  be  small  so  that  photoelectrons  are  released  from  a 
small  area  on  the  detector.  This  physical  reasoning  is 
described  by  the  Y parameter. 

With  regard  to  the  above  discussion  and  by  careful 
observation  of  Eas  (26)  and  (27),  the  upper  and  lower  bounds 
converge  to  the  steady  state  value  of  the  actual  MSE  when 
X>>1  and  Y<<1.  Referring  back  to  the  performance  parameters, 
the  estimator  will  perform  well  if  a large  number  of  signal 
photoelectrons  are  observed  frequently  and  they  appear 
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in  a small  area  of  a photodetector  surface.  This  physical 
reasoning  is  described  by  the  above  parameters  derived  from 


upper  and  lower  bounds.  These  parameters  establish  regions  of 
operation  where  the  upper  and  lower  bounds  converge  to  the  actual 
mean  square  error.  However,  when  the  bounds  do  not  converge, 
actual  performance  is  difficult  to  obtain.  To  demonstrate 
the  usefulness  of  the  bounds  and  to  determine  which  bound  lies 
closer  to  actual  performanc a Monte  Carlo  simulation  is 
conducted.  These  parameters  and  the  Snyder  Filter  assumes 
no  dark  current  in  the  measurement  process;  and  so  an  additional 
objective  is  to  include  d.ark  current  in  the  simulation  to 
determine  its  effects  in  filter  performance. 

Summary 

The  chapter  began  by  describing  the  direct  detection  optical 
receiver.  The  important  point  to  remember  is  that  the  statistics 
of  the  output  of  the  receiver  govern  the  performance  of  optical 
components  performing  postdetection  processing  and  vice  versa. 

These  statistics  are  described  by  a mathematical  model  called 
the  point  process.  The  model  is  applied  in  situations  when  the 
receiver  is  operating  under  low  signal  and  background  radiation. 

In  this  case,  the  discrete  nature  of  light  at  optical  frequencies 
and  the  nature  of  the  detection  process  fundamentally  limit  the 
performance  of  optical  trackers. 

The  Snyder  Filter  was  discussed  next.  This  filter,  a recursive 
data  processing  algorithm,  processes  the  point  process  observations 
of  the  receiver  to  estimate  optimally  the  position  of 
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maximum  intensity  on  the  detector  surface.  The  Snyder  Filter 
assumes  its  observations  are  from  an  infinitesimally  fine 
detector  array  (i.e.,  the  detector  array  has  infinite  resolution 
to  measure  the  exact  location  of  a photoconversion). 

The  chapter  concluded  with  the  performance  parameters 
describing  the  mean-square-error  of  the  Snyder  Filter.  The 
parameters  are  derived  from  the  upper  and  lower  bound  on  MSE. 

The  two  parameters  are  defined  to  establish  conditions  under 
which  the  Snyder  Filter  should  yield  good  tracking.  One  parameter 
is  the  ensemble  average  of  the  number  of  photons  detected  in  a 
coherence  time  of  the  spot's  dynamics,  and  the  other  is  the 
spot's  jitter  normalized  by  the  square  of  the  beam  width.  These 
parameters  also  indicate  sufficient  conditions,  such  that  the 
bounds  converge  to  the  actual  MSE.  The  MSE  can  then  be  used 
to  evaluate  the  tracking  capabilities  of  the  Snyder  Filter 
(i.e.,  to  determine  whether  the  average  errors  committed 
by  the  Snyder  Filter  are  sufficiently  small  for  tracking  purposes). 
Again,  these  performance  parameters  and  the  Snyder  Filter 
assume  that  the  events  locations  are  only  due  to  the  spot. 
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III.  Simulation  of  the  Snyder  Filter  and  Performance  Analysis 

The  last  chapter  addressed  the  following  concepts:  (1) 
a mathematical  model  of  the  output  process  of  the  optical 
receiver,  (2)  a model  of  the  spot's  intensity  profile  and 
its  dynamics,  (3)  the  Snyder  Filter  which  incorporates  the 
measurement  data  { t^ } and  {r^}  from  the  output  of  the 
receiver  to  estimate  the  position  of  maximum  intensity,  and 
(4)  performance  parameters  which  establish  regions  of  opera- 
tions where  the  bounds  converge  to  actual  MSE. 

This  chapter  develops  algorithms  to  implement  the 
above  models  and  to  simulate  the  Snyder  Filter  in  a digital 
computer.  The  algorithms  that  were  developed  are  (1) 
simulation  of  spot's  dynamics,  (2)  a noise  generator,  and 
(3)  algorithms  to  generate  the  time-space  point  process 
observation  (two  methods). 

Concluding  the  chapter  is  the  performance  analysis  used 
to  evaluate  the  tracking  capabilities  of  the  Snyder  Filter. 

Simulation  of  Spot's  Dynamics 

The  spot’s  dynamics  on  the  detector  surface  modelled 
by  a linear  stochastic  differential  equation  given  in  Eq.  (12), 
or  in  less  rigorous  form, 

x(t)  = F(t)x(t)  + G(t)w(t)  (30) 

where  w(t)  is  a white  noise  with  zero  mean  and  variance 
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E[w(t)w(t»  )]  = Q ( t ) <5  (t-t ' ) 


where  Q(t)  is  a diffusion  strength  variable.  Consider  an 
equivalent  discrete-time  model  for  digital  simulation.  Here, 
equivalence  means  that  a discrete  time  process  has  the  same 
second  order  statistics  as  x(t),  when  viewed  at  the  sample 
times.  As  shown  in  Ref.  (9: 4.63),  an  equivalent  discrete 
time  solution  to  Eq.  (30)  is 

x(t1+i)  = + ( t i ) (32) 

where  t1+^,td  are  the  discrete  points  in  time,  4>  ( t i+  * t i ) is  the 
state  transition  matrix  (scalar,  in  this  case),  w^(t^)  is  the 
white  discrete-time  stochastic  process  with  statistics: 


E[wd(ti)]  = 0 


E[wd2(t1)]  = Qd(t1) 


t . 

= / + 1<j>  (t1+1,  T )G  (t  )Q(  t )G  ( T )<|>  (t  +1,  T )dt 

h • t = t 

, tj 


E[w(J(t1)wd(tj  )]  = 0 


: h * ‘j 


Prot  Ref.  (9;30),the  mean  m (t)  and  variance  P (t)  of  the 

' X XX 

process  x(t)  propagate  in  time  as 
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Pxx^ti+l)  ~ <t(ti+l,ti^Pxx(ti)4>(ti+l»ti)  + Qd(ti^  (3i|) 

A simple  model  of  the  spot's  dynamics  on  the  detector  is 
an  output  process  of  a first  order  lag  driven  by  white 
Gaussian  noise  of  zero  mean.  This  model  is  represented  as 

x(t)  = - -i-  x ( t ) + w(t)  (35) 

Tc 

where  F(t)  = - — and  G(t)  = 1.  Then,  its  state  transition 
Tc 

matrix  is 

4>  ( t , t ) = <J>  ( t — t ) = e tc^ 

The  model  therefore  becomes. 

At 

x(ti+i)  = e tc  xCtj^)  + wd ( 1 1 ) 

where  At  is  sample  period  (or  observation  time  increment), 
and  wd(fcd)  is  a white  Gaussian  discrete-time  process  of  zero 
mean  and  covariance 


, 2At 

= •/  ti  ■ 


2 2 2 

where  Q = — o , and  o 4 mean  square  value  of  output 

T _ X X = 

c 2 

process  x(t),  (in  this  report,  a , represents  the  "jitter" 
of  spot).  Q is  found  in  Ref.  9:^-77  ; and  it  represents  the 
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the  strength  of  the  white  noise.  Eq.  (35)  is  the  expression 
which  simulates  the  spot's  dynamics  on  the  detector.  The 
next  section  discusses  how  a white  Gaussian  discrete-time 
process  w^(t)  is  generated. 

Noise  Generators 

A uniform  random  variable  distributed  between  0 and  1 
is  supplied  by  the  computer  system.  This  uniform  number 
generator  is  used  to  generate  realizations  of  random  variables 
described  by  desired  probability  density  functions  or  proba- 
bility distribution  functions.  This  will  be  shown  in  the 
next  section. 

However,  applying  the  Central  Limit  of  Statistics,  one 
can  generate  an  essentially  Gaussian  random  variable  by 
adding  twelve  independent  uniform  random  variables  distributed 
between  0 and  1.  Since  the  mean  and  variance  of  a single 
uniform  random  variable  is  1/2  and  1/12,  respectively,  the 
above  addition  yields  a Gaussian  random  variable  with  mean  6 
and  variance  1.  To  obtain  a zero  mean  and  a unit  variance 
simply  subtract  6 from  the  above  result.  Note  that  the 
Central  Limit  Theorem  implies  that  the  resulting  random 
variable  w has  a probability  distribution  that  approximates 
a Gaussian  distribution  (Ref.  11). 

Finally,  a Gaussian  random  variable  w.  of  desired 
variance  ow  and  mean  v>w  is  given  by  the  following  transfor- 
mation equation 
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w = o z + y 
w w 

where  z is  a Gaussian  random  variable  with  zero  mean  and  unit 
variance . 

The  "whiteness"  of  the  Gaussian  discrete-time  sequence 
wd(t^),  is  obtained  by  independent  calls  of  the  above  algorithm 
The  resulting  set  of  w's  forms  the  white-Gauss ian  discrete- 
time process  w(J(t1)  (i.e.,  E [ ( t ^ ) w^ ( t j ) ] = 0 for  all 
tf  5*  tj).  Also,  note  that  the  second  order  statistics  of 
wd(tf)  (or  set  of  random  variables  w)  are  the  second  order 
statistics  of  the  output  process  x(t),  of  the  first  order 
lag.  The  statistics  of  w^(t^)  are 


? 2At 
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Point  Process  Generation  in  Time  and  Space 

In  this  report,  two  methods  are  employed  to  generate 
the  event  locations  in  time  and  space.  One,  the  Binomial 
Method,  is  based  on  the  incremental  definition  of  the  point 
process,  and  the  other,  the  Conditional  Method,  makes  use 
of  conditional  probability  densities  describing  the  output 
statistics  of  the  receiver.  The  latter  method  uses  Bayes* 
Rule.  It  should  be  noted  that  these  conditional  densities 
are  derived  by  the  incremental  definition  of  a point  process. 
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The  Binomial  rule  will  be  discussed  first.  Next,  a trans- 
formation equation  used  by  the  Conditional  Method  is  described; 
then  the  Conditional  Method  follows.  Concluding  this  section 
is  a discussion  on  the  theoretical  and  practical  aspects  of 
the  two  methods  of  simulation. 

Binomial  Method . Shown  in  Pig.  (7)  is  a flowchart  of  the 
Binomial  algorithm  that  generates  the  time-space  events.  As 
mentioned  earlier,  this  development  is  based  on  the  incremental 
definition  of  a point  process.  It  is  the  same  development 
describing  the  model  of  the  photodetecting  surface.  For 
the  purpose  of  clarity,  a discussion  on  the  point  process 
incremental  definition  in  one-dimensional  space  is  given  as 

(1)  Event  in  a space-time  observation  area,  AX. At,  are 
statistically  independent  where  AX,  is  the  incremental  detector 
length.  At  is  the  incremental  time  observation. 


P 


one  event  in 
A£  At 

located  at  (t,r) 


= X(r,x(t)  )AtAX,  + o(  At  AX. ) 


(36) 


P 


zero  events  in 
AXAt 

located  at  (t,r) 


1 - X(r,x(t) )AtAX  + o( At  AX ) (37) 


where  P[»]  denotes  probability  of,  the  " oAXAt"  denotes  higher 
order  terms  of  AXAt  (i.e.,  ( ° [ ^ ) = 0),  and  x(t)  denotes 

location  of  maximum  intensity  on  the  detector  surface.  Note 
that  different  values  of  r will  give  different  values  of  rate 
functions,  X(r,x(t). 

Consider  a uniform  random  number  distributed  between  zero  and 
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Figure  7.  Flowchart  of  Binomial  Method 


one.  Then,  the  probability  of  a value  between  zero  and  a 

number  \p , is  simply  ip . So  let  ip  = X ( r,x(t ) ) ) At.At  establish 

a threshold  to  declare  an  event, where  U„  is  called  from  a 

k 

uniform  random  generator  (that  is,u  < X (r , x ( t ) ) At /u  means  an 

k 

event  occurred  at  r on  the  detector  surface). 

The  above  module  takes  into  account  a different  threshold 
and  rate  function  when  moving  from  one  spatial  point  to 
another  on  the  detector  length.  This  method  also  includes 
effects  of  dark  current  (i.e.,  X(r,x(t))  = X (r,x(t))  + X^, 


where  A^4  Gaussian  rate  function  \ 4 dark  rate  function).  The 
Binomial  Method  imposes  a significant  computational  load  in  the 
simulation,  as  will  be  shown  in  a later  section  following  the 
Conditional  Method.  Note  that  this  method  is  an  approximation 
for  finite  M and  At  and  that  there  is  an  issue  of  coarseness 
of  the  time  and  space  discretization.  This  creates  implemen- 
tation and  theoretical  aspects  which  will  be  discussed  in  a 
later  section. 

Transformation  Equaticn.  The  transformation  equation 
is  used  to  generate  a desired  random  variable  of  a given 
probability  distribution  function  by  calling  a uniform  random 
generator.  In  this  case,  either  an  event  location  or  event 
time  is  generated  by  calling  a uniform  random  variable  U 
which  is  distributed  between  zero  and  one  with  probability 
density  height  1.  The  solution  to  this  transformation  problem 
is  found  in  Ref.  (12:184)  and  is  given  as 

U = Fd(D)  (38) 

or 

D = ?D-1(U)  (39) 

where  F^  is  the  desired  distribution  function,  D is  the  desired 
random  variable,  F^-^  is  the  inverse  function  of  F^ . Fig.  (8) 
graphically  depicts  this  solution. 

The  random  variable.  D represents  either  the  random  event 
location  or  time  of  the  time-space  point  process.  Similarly, 

Fd  (or  equivalently,  the  transformation  function),  corresponds 
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Figure  8.  Generation  of  a Specified  Probability 

Distribution  Function  (Copied  from  Ref.  12:184) 

to  either  the  event  location  or  time  probability  distribution. 

In  this  report,  these  probability  distributions  are  derived  from 
the  densities  f[rk|t^,n]  and  f[sk]  where  s^  is  the  interarrival 
times  of  ordered  events  1,  (k=0,  l,2...,n),  t^  is  the  event 
times  from  a time  orgin,  and  n is  the  number  of  events.  These 
densities  are  shown  in  the  next  subsection. 

Conditional  Method . To  generate  the  time-space  point 
process  events  from  the  output  process  of  the  detector  surface, 


the  number  of  events,  the  event  times  and  event  locations 
are  needed.  Here,  the  probability  density  f [ { rk} , { t^ } , n]  is  required 
to  characterize  output. This  joint  density  is  employed,  to  derive 
the  conditional  densities  of  event  times  and  event  locations 
conditioned  on  knowing  certain  variables, where  the  variables 
are  generated  by  actual  implementation  of  the  simulation.  The 
development  of  this  section  will  involve  some  lengthy  (but  not 
complex)  derivations  of  expressions  needed  to  generate  {r^.}  and 
{t^}.  Also,  physical  insights  about  the  Snyder  Filter's  tracking 
capabilities  will  be  cited  where  appropriate. 

In  this  study,  the  variables  IQ(t),a(t),  F(t),  G(t), 

H(t),  and  Q(t)  are  time-invariant.  Then  the  average  rate  at 
which  photoelectrons  appear  on  the  photodetectors  array  (after 
integrating  the  spatial  dependence  of  the  rate  function, 
defined  in  Eq.  (8),  over  entire  detector  length)  yields 

X0  = *0^  0 + Xd£ 

where  Eq.  (14)  was  used,  and  l is  the  detector  length  of  array. 
Eq.(40)  has  terms  including  the  effects  due  to  dark  current; 
however,  in  this  relation,  some  assumptions  were  used,  brought 
forth  earlier  in  the  report.  To  explain  the  first  term  in  Eq . 

(40),  it  is  assumed  that  the  detector  surface  is  large,  such 
that,  with  respect  to  the  spot's  Gaussian  intensity  profile, 
the  detector  surface  is  infintely  large.  Thus,  edge  effects  of 
the  spot's  intensity  profile  on  the  detector  may  be  Ignored. 

On  the  other  hand,  \ is  finite  which  means  that  /“  X.dj, 

d loo  d 
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_2 

2 

has  an  infinite  value.  Therefore,  the  integral,  /_*  xdd«. 

2 

results  in  the  second  term  of  Eq.  (40).  The  major  point  is 
that  the  detector  length  must  be  large  so  that  edge  effects 
between  the  spot  and  detector  are  neglible  when  the  spot  is 
jittering  about  on  the  detector  surface. 

Consider  the  simulation  of  event  times.  It  is  known 
that  the  interarrival  time,  denoted  as  s^,  has  an  exponential 
density  function  (Ref.  3:60)  given  as 

f[sk]  = X0exp(-XQs)  (4l) 

To  obtain  integrate  Eq.  (4l)  over  the  time  continuum 

between  and  s^.  Using  the  transformation  function  from 
Eq.  (38),  Eq.  (4l)  yields 

V = F[sk]  = 1 “ exp(-A0sk)  (42) 

where  U^'  is  a uniform  random  variable  distributed  between 
zero  and  one.  It  then  follows  that  Uk  = 1 - U^'  is  also  a 
uniform  variable  with  the  same  statistics.  By  straight  forward 
manipulation,  Eq.  (42)  yields 

sk  ' -7^  tn  \ ("3) 

To  get  the  event  times,  just  use  a running  sum.  In  this  case, 
with  tg  as  the  initial  time  the  event  times,  of  the  point  process 
with  rate  function  Xq  can  then  be  simulated  by  the  following 
algorithm:  call  a uniform  random  variable,  use  Eq.  (43)  to 

transform  to  s^;  assign  s^  as  the  time  from  tg  to  the  first 
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photoconversion  (i.e,,  = tg+s^);  call  another  uniform 

random  variable,  Up;  transform  Up  to  s2>  assign  Sp  as  the  time 
between  the  first  and  second  photoconversion  (i.e.,  Tp  = + 

s?),  and  so  forth  (Ref.  3:62). 

The  next  problem  is  finding  an  algorithm  to  generate 
the  event  locations,  r^'s.  This  problem  is  solved  by  using 
densities  f [ {r^} , { t^ } , n]  and  f[{t,^.},n]  (see  Appendix  A for 
derivation).  These  densities  describe  the  statistics  of  the 
output  process  from  the  photodetector  surface  which  are  expressed 
as 

f[(t.  },  (r.  },n}]  = exp ( J {T  X (a , & )dadB  ) II  X(t.  ,r,  ) (44) 

K Ad  u k=l  K k 


f[{tk),n]  = exp ( - Jq  X ( a ) da ) n X C t^ ) 

k=l 


(45) 


where 


X(tk,rk)  = -f-  exp (- 


nl«  -(r,-x(t.))‘ 


) + X, 


2 a 


VW 


n I 


A(tk^  ~ fQ  ^ 0 + XD 


(46) 


(47) 


Note  that  f[{tk>,n]  is  the  marginal  density  of  f [ { t^ } , {r^ } , n] 
and  that  X(t^)  in  Eq.  (47)  has  no  spatial  or  temporal  dependence. 
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Also,  the  integrals  in  Eq.  (44)  and  (45)  will  both  yield 
AqT . Thus  Eqs . (44)  and  (45)  respectively  yield 


f[(tk), (rk),n]  = exp ( -XqT)  n l(tk,rk) 

k ~ 1 


f[{t,  },n]  = exp(-XnT)  n X(t,  ) 
K U k=l  K 


With  the  help  of  Bayes'  rule,  the  conditional  probability 
density  function  of  event  locations  conditioned  on  knowing  the 
event  times  and  the  number  of  events  is  given  as 


f[{rk}  I {tk}’n]  = 


f[{tk>,n] 


n WV  + xd 


= n 
k=l 


f[{i’k)  | {tk},n]  = 


X (t^,r^)  + X^  X 2,r  2^  + ^ 


n >gCtk.rk)  ^ >d 

1=3  x0 


fl^ri  I tj_,n]  f 2 I 


n,  fk[r'k|tk>n:| 


58 


where  Eqs.  (48)  and  (49)  were  used.  The  above  relation  says 
that  the  joint  conditional  density  f [ {r^ } | { t^} ,n] , can 
be  derived  by  forming  the  product  of  each  of  the  densities  at 
the  time-space  coordinates,  C t ^ ) • Thus,  for  n=l  at  (r^t^), 
the  conditional  density  of  r^  given  t^  is  expressed  as 


f[rk|tk]  - kiVV 

Ao 


Before  proceeding  further,  it  is  useful  to  check  if  Eq.  (51) 
is  normalized  (i.e.,  f [a^  | t^Jda^l)  . That  is. 


-(a,  -x) 

-fr  Jo  exp  2 M<*k 


= _L  (_H_ 

xo  *lfo 


Iq/27  a 


+ xdO 


1 


Again,  it  is  assumed  that  the  detector  is  large  such  that  the 

limits  in  the  integral  of  the  first  term  on  the  left  side  of 

Eq.(52)  are  approximately  valid. 

To  get  the  probability  distribution  function  F[r^|t^] 

integrate  f[rk|t^]  over  the  detector  plane  from  -°°  to  r^  for 

£ 

the  first  term  and  from  — ^ to  r^  for  the  second  term  in  Eq . 
(51).  The  result  is  expressed  as 


erf  v A j — exp(-Z  jdZ 

* 7T  ' (J 


rk-x(tk) 

r r « 


Using  the  transformation  equation  defined  in  Eq.  (39)  <ind  after 
algebraic  manipulation,  Eq.  (53)  yields 


erf  v + (—  Ap  v = 2A2  Uk  - 1 + 2A,,  x(tfc)  + | (54) 

/rr 

where 


A,  A — =r 

1 = n I, 


A2  = (•'  fQ  /2v  o + X l) 


hfg  /2tt  o 


A3  = 


£ /27  a 


Eq.  (54)  is  used  to  generate  the  event  locations  at  the 
event  times  given  the  variables  defined  in  this  expression, 
The  event  locations  are  generated  by  a table  lookup  method 
with  linear  interpolation  and  is  performed  as  follows: 

(1)  Evaluate  and  obtain  all  necessary  variables, 
i.e.,  A-p  A2,  A^,  -j,  U^,  and  x(t^). 
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(2)  For  a specific  range  of  values  for  v,  solve  the 
expression  on  the  left-hand  side  of  Eq.  (54),  and  store  the 
results  in  a table  (i.e.,  computer  memory ), 

(3)  Derive  a random  variable  s by  calling  a 
uniform  random  variable,  Uk,  distributed  between  (0,1)  and 
by  obtaining  the  true  position  of  the  centroid,  x(t),  from 
the  truth  model  simulation.  Given  these  values,  the  right- 
hand  side  of  Eq.  (54)  is  evaluated  and  the  random  variable 
realization  s is  obtained, 

(4)  To  find  v,  obtain  s by  a table  lookup  with  linear 
interpolation.  In  the  table  the  increment  of  v is  equal  to 
0.025,  where  erf  v is  piecewise  linear  for  this  approximation; 
thus,  linear  interpolation  is  used 

(5)  Upon  finding  v,  it  then  follows  that  the  event  location 
is  evaluated  by  the  following  expression: 

rk  = (/"2~  o)v  + x(tk) 

This  paragraph  completes  the  algorithm  of  generating  the 
event  locations. 

Since  Eq.  (54)  governs  the  generation  of  event  locations, 
it  is  useful  to  examine  this  relation  for  limiting  cases 
(i.e.,  1^=0  and  X (t;k,rk)  = 0).  This  investigation  will 
provide  insights  into  the  tracking  capabilities  of  the  Snyder 
Filter . 

= 0 implies  A^  = A^  = 0 and  A2  = 1.  Eq.  (54) 
is  then  manipulated  into  a convenient  form  as 
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(55) 


r 


/ /rk_X(tk)XN 

U = l/2(l  + erf  (-£ £-)) 

K /2  a 


uk  ■ VhJV 


(56) 


where  Eq.  (38)  is  implemented  in  tjie  last  expression.  Corresponding 
to  (r^.  1 1^.)  is  f(rjc|t^)  which  is  given  as 

f[rk|tk]  = — — exp{-^-[rk-x(tk)]2}2  (57) 

/2ttct  2 a 

The  event  locations  are  generated  by  calling  a Gaussian 
random  variable  of  mean  x(tk)  and  variance  a,  at  the  event 
times.  Note  that  the  intensity  profile  of  the  spot  being 
Gaussian-shaped  leads  to  a result  that  the  event  locations  are 
generated  as  Gaussian  random  variables  conditioned  on  knowing 
the  event  times.  The  Snyder  Filter  is  based  on  the  assumption 
of  a Gaussian-shaped  profile  of  the  spot.  Thus,  for  A^  / 0, 
"smears"  out  the  Gaussian  information  as  seen  in  Eq.  (54). 

For  Ag (r , t ) >> A^ , the  effect  of  dark  current  is  examined 
by  doing  a perturbation  analysis.  This  analysis  is  performed 
by  factoring  out  the  larger  of  the  two  rate  functions  in  Eq . 

(47).  Then  the  natural  logarithm  of  the  left  and  right  hand 
side  of  Eq . (47)  becomes 

n A,  n A (r,  ,t,  ) 
in  f({rk>|{tk>,n)  » £ *n-r^  + £ — lL_ 

K k=l  A0  k=l  A0 


/ ^g^rk’^k^\  VVV 

where  the  approximation  In  ( 1 + ^ ) = -■»—  - — — 

V Ad  / Ad 


and 
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<<  1 Is  used.  The  exponential  of  the  left  and 


V 


right  hand  side  of  the  above  equation  yields 

exp  ( 


f[{r>k)  |{tk),n]  a x 


n 

n 


»d  A(rk>V 


(58) 


Again,  Eq  (58),  can  be  derived  by  forming  the  product  of  the 
above  individual  densities  at  (rk,tk).  Then,  setting  Xg(rk,tk)  = 0 
in  Equation  (58)  and  applying  this  condition  to  Eq  (40), 

Eq  (58)  yields 


f[rkitk]  - 


1 

£ 


(59) 


where  n = 1.  Intuitively,  Eq  (58)  is  reasonable,  since  at 

any  given  time  a dark  event  can  occur  in  any  differential 

length  on  the  detector  surface  with  equal  probability. 

In  summary,  for  signal-limited  conditions,  the  event  locations 

are  called  as  a Gaussian  random  variable  with  mean  x(tk) 

2 

and  variance  o ; for  dark-limited  conditions  the  event  locations 
are  generated  as  uniform  random  variables  distributed  in 
(-  2>2^  with  probability  density  height  — ; and  for  generating 
event  locations  other  than  the  limiting  cases,  Eq  (54)  is 
implemented.  Finally,  the  generation  of  event  locations  is 
called  at  the  event  times  given  by  Eq  (43). 

Aspects  of  Implementing  each  method . The  Binomial  Method 
has  two  restrictions  for  its  algorithm  to  work  properly: 
one,  is  x[r  ,x(t  )]A£At<l  and  the  other,  is  A£At<<l  where  for 
implementation  purposes  At  4 sample  time  or  observation  time 


63 


interval.  Because  of  these  conditions,  several  practical  and 
theoretical  problems  arise.  One  is  a resolution  problem.  The 
Snyder  Filter  assumes  that  its  position-sensitive  measurements  are 
from  an  infinitesimally  fire  partitioned  detector  array.  No 
matter  how  small  " M"  is,  there  is  still  a problem  in 
extracting  the  "exact"  location  of  the  event.  For  this  report, 
the  detector  array  is  divided  into  equal  increments  of  A l,  and 
the  value  of  the  rate  function  at  the  center  of  each  incremental 
length,  is  chosen  to  evaluate  the  probability,  X (r ,x(t ) ) AHAt . 

This  procedure  determines  where  an  event  (if  any)  occurred 
on  the  detector  surface.  This  method  is  depicted  in  Fig  (9). 

Note  that  At  is  the  observation  increment  during  the  simulation 
time.  From  a practical  point  of  view,  this  method  consumes 
a large  amount  of  computer  time; at  each  time  increment  of 
observation ,A£ is  scanned  across  the  detector  surface  from  one 
end  to  the  other. 

In  the  Binomial  Method,  a theoretical  point  that  should 
be  considered  is  the  question  of  "simultaneous"  events  occurring 
on  the  surface.  The  Snyder  Filter  assumes  that  one  event  is 
observed  at  any  one  time.  With  this  simulation  method,,  there 
is  a probability  that  more  than  one  event  can  occur.  Even 
ifAt<<l  and  A a< <1 , there  is  still  a probability  (though  small) 
that  simultaneous  events  can  occur  at  different  locations. 

It  was  mentioned  previously  that  scanning  the  detector  surface 
for  small  a£  results  in  a large  comsumption  of  computer  time 
at  each  time  increment.  Thus,  if  a smaller  At  or  AJl  is  chosen 
to  reduce  the  probability  of  simultaneous  events,  this  will 
result  in  a larger  consumption  of  computer  resources. 
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Fig.  9.  Resolution  problem  in  Binomial  Method 

Furthermore,  when  dark  current  is  included  in  the  simulation, 
the  probability  of  simultaneous  events  grows  substantially. 

To  reduce  the  amount  of  computer  resources  (and  to  observe  the 
errors  generated  by  the  Snyder  Filter  between  event  times), 
the  rate  functions  are  chosen  to  have  a low  value.  This 
enables  the  use  of  a larger  sample  period.  At,  which  also 
implies  a higher  probability  of  simultaneous  events.  It  is 
evident  now  that  the  Binomial  method  has  implementation 
problems.  In  spite  of  these  problems  this  method  is  performed 
in  the  computer,  and  the  results  are  shown  in  the  next  chapter. 

In  terms  of  implementation,  the  Conditional  Method  is  a 
better  algorithm,  since  the  event  locations  are  called  at  the 
event  times;  (i.e.,  the  problem  of  simultaneous  events  is 
precluded).  In  the  simulation  it  is  desirable  to  sample 
the. errors  committed  by  the  Snyder  Filter  between  successive 
event  times.  In  this  case,  it  is  convenient  to  break  up  the 
simulation  time  interval  (observation  time)  into  time  increments. 
At.  At  will  denote  the  sample  period  in  order  to  sample  the 
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error  process  (as  done  in  the  Binomial  Method).  When  the  event 
times  are  generated,  there  is  a situation  that  the  interarrival 


times,  w^'s,  can  occur  smaller  than  the  sample  period  (i.e. 
s^  < At).  This  situation  is  shown  in  figure  (10). 

Performance  Analysis 

The  objective  of  the  performance  analysis  is  to  characterize 
the  error  process  statistically.  This  is  done  by  generating 
many  samples  of  the  error  process  through  Monte-Carlo  simulation. 
The  simulation  uses  random  number  generators  and  shaping  filters 
to  generate  the  random  errors.  One  objective  in  the  simulation 
is  to  portray  filter  performance  capabilities  when  parameters 
that  describe  the  problem  are  varied.  First,  consider  the  test 
cases  without  dark  current.  Here,  the  signal  rate  function  or 
spot's  jitter  are  changed  and, thus,  filter  capabilities  can 
be  evaluated.  Note  that  changing  the  signal  rate  varies  the 
average  number  of  photons  detected  per  coherence  time;  similarly, 
changing  the  spot's  jitter  causes  variations  in  the  Y parameter. 

In  this  study,  the  test  cases  do  not  meet  sufficient  conditions 
which  are  X>>1  and  Y<<1  ; therefore, the  results  should  provide 
an  indication  of  which  bound  is  close  to  actual  performance. 

A plot  which  contains  the  sample  variance  and  average  filter 
variance  is  a useful  output  to  make  this  comparison. 

Second,  the  measurement  model  (of  wh  ich  the  Snyder  Filter 
is  based)  is  altered  by  presenting  dark  events  in  the  measurement 
process.  The  above  plot  is  again  a useful  output  to  evaluate 
the  filter's  sensitivity  to  dark  measurements. 

In  this  study,  a final  objective  of  the  performance  analysis 
is  to  change  the  filter  structure,  through  ad  hoc  methods,  when 
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Fig  10.  Sampling  the  Error  Process  with  Fixed  Sample  Period 
(Conditional  Method) 


satisfactory  performance  cannot  be  obtained  in  the  presence  of 
dark  events.  Again,  the  above  plot  containing  the  variances 
would  be  useful  in  providing  insights  of  how  to  alter  the- 
ft Iter  structure. 

As  previously  mentioned,  the  error  process  is  characterized  by 
generating  many  sample  runs  of  this  process.  The  appropriate 
number  should  be  chosen  sufficiently  large,  such  that  additional 
sample  runs  does  not  appreciably  change  the  sample  statistics. 

The  last  statement  is  an  illustration  of  the  weak  law  of  large 
numbers:  let  E[e]  = e denote  the  ensemble  average  of  the 

error  process  e.  Also,  let  <e>^  denote  the  empirical  average  for 
N runs.  When  N is  large,  the  weak  law  states  with  high 


probability  that  <e>..  is  close  to  the  number  e.  In  mathematical  terms, 

a-2 

r[|<e>„  - e|>E]  = — 


- 1 

where  <e>„  4 ^ J e.  sample  mean  error,  N 4 number  of 
N j = l J 

statistical  independent  random  variables,  e ^ with  mean,  e 

2 

and  variance  0—  , and  e is  an  arbitrary  constant.  The 


weak  law  Is  a measure  of  the  confidence  of  inferriny  Information 
about  the  actual  statistics  of  the  error  process  (Ref  11:96). 
Thus,  as  the  number  of  sample  runs  increases,  the  confidence 
in  characterizing  the  error  process  accurately  increases. 


Sample  statistics . The  statistics  computed  at  each  time 

point  are  the  calculated  mean  error  e,  calculated  variance 
2 

CTe  , and  the  ensemble  average  of  the  filter  variance  terms 
PF.  These  calculations  are  performed  over  the  ensemble  of 
runs  N,  for  each  time  point  tm<  These  expressions  are 


N 


M 1 [X(t  ) 

N L L m 

J-l 


x ( t ) ] 
m 


N 

P £[x  t ) -x(t  )]‘ 

2ctj  - -4-  m 


'§  ' m'  j = l N-l 


N N[<e(t  )>„]' 
£ m N 


J-l  N-l 


N 


PF  (t)  = h 1 PF ( t ) 
m Nj“i  m 


(60) 


(61) 


(62) 


The  above  relations  are  found  in  Ref  (13).  Note  that  the  time 
point  t are  not  event  times  but  are  the  sample  times 
for  a fixed  sample  period  At.  From  Eq  (6l)f  the  variance 

p 

Og.  (t  ) is  divided  by  (N-l)  instead  of  N?  which  results  in  an 
unbiased  estimator  of  variance  (Ref  13). 

Plots . Mentioned  in  previous  paragraphs  is  that  a useful 
output  in  the  analysis  is  a plot  of  the  average  filter  estimate 
of  the  error  variance,  along  with  the  corresponding  computed 
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rr  ^ 

sample  Cor  true)  variance  for  all  times  of  interest.  This  output 
is  used  to  compare  if  the  sample  variance  of  the  error  process 
matches  with  the  filter's  average  error  variance.  In  other 
words,  this  output  should  indicate  how  well  the  filter  is 
estimating  its  errors  with  the  actual  errors  it  commits;  and 
to  compare  with  the  bounds  on  MSE  for  test  cases  where  dark 
current  is  set  equal  to  zero. 

A second  useful  output  on  this  analysis  is  a plot  of  the 
mean  error  versus  time.  Along  with  this  plot  are  the  standard 
deviation  bounds  plotted  above  and  below  the  mean  error.  This 
plot  is  useful  if  the  Snyder  Filter  provides  the  desired  accuracy. 

The  two  plots  that  were  just  described,  are  used  to  analyze  the 
average  error  performance  of  the  Snyder  Filter.  For  all 
cases,  these  plots  are  used  to  determine  whether  the  actual 
errors  committed  by  the  filter  are  sufficiently  small  to  yield 
good  tracking  performance. 

Now  consider  plots  of  individual  sample  of  the  error 
process.  For  no  dark  current,  this  plot  is  used  to  observe  how 
the  Snyder  Filter  minimizes  its  errors.  When  dark  current  is 
represented  in  the  simulation,  then  the  Snyder  Filter  sensitivity 
to  the  dark  measurements  can  be  observed.  Along  the  zero  axis 
of  this  plot  are  the  event  times  to  indicate  clearly  how  the 
filter  updates  its  estimate  at  the  event  times.  Also,  a plot 
of  the  sample  filter  variance  is  also  of  interest,  in  order  to 
observe  how  the  Snyder  Filter  propagates  the  uncertainty  of  Its 
estimate  in  time. 
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Test  cases. 


During  this  study,  different  cases  using 


i 


different  values  of  X and  Y performance  parameters  are  considered. 
The  test  cases  are  given  in  tables  (1)  and  (2).  All  test  cases 
are  performed  where  there  is  a relative  separation  between  the 
bounds.  Here^  the  bound  lying  closer  to  actual  performance  can 
be  determined.  For  the  test  cases  in  Table  (1)  and  (2),  the 
dark  current  is  set  equal  to  zero. 

When  dark  current  is  presented  in  the  simulation  separate 
test  cases  are  performed.  The  measurement  process  of  the  Snyder 
Filter  is  altered  by  this  addition  of  dark  current.  For  these 
test  cases,  a sensitivity  parameter  is  defined  to  be  the  average 
signal-to-noise  counts  ratio  which  is  expressed  as 


S 

N 


Here,  the  degree  of  degradation  in  filter  performance  is 
observed  as  the  average  signal-to-noise  counts  ratio  is 
varied.  Because  of  the  time  constraints  imposed  in  the  study, 
limited  cases  for  Xd  ^ 0 are  examined. 


Summary 

The  chapter  developed  the  necessary  algorithms  needed  to 
perform  this  study.  These  algorithms  consists  of  simulating  the 
spot's  dynamics,  and  the  generation  of  position-sensitive 
measurements  (i.e.  time-space  point  process).  The  chapter 
concluded  with  a general  discussion  on  performance  analysis, 
where  it  consists  of  sample  statistics,  plots,  and  test  cases 
needed  to  investigate  the  tracking  capabilities  of  the  Snyder  Filter. 
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Performance  Parameter 


X 

Y 

5.0 

0.1 

5.0 

0.25 

5.0 

1.0 

Table  1.  Test  Cases  for  Binomial  Method 
(without  dark  current) 


Performance  Parameter 


X 

Y 

5.0 

0.5 

5.0 

1.0 

5.0 

5.0 

5.0 

10.0 

1.0 

0.5 

1.0 

1.0 

1.0 

5.0 

1.0 

10.0 

0.5 

5.0 

10.0 

5.0 

15.0 

5.0 

20.0 

5.0 

0.5 

1.0 

10.0 

1.0 

15.0 

1.0 

20.0 

1.0 

Table  2.  Test  Cases  for  Conditional  Method 
(without  dark  current) 
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IV  Results 


; 

This  chapter  contains  results  from  the  Monte  Carlo  simulation. 

The  simulation  will  indicate  how  well  the  actual  performance  can 
be  described,  when  compared  to  the  upper  and  lower  bounds.  Thus, 
various  values  on  the  set  of  parameters,  derived  from  the  bounds, 
are  used  to  determine  which  bound  is  closer  to  actual  performance. 
The  bounds  are  valid  only  when  the  measurement  is  perfect  (i.e., 
there  are  no  dark  or  background  noise  events).  When  dark  current 
is  included  in  the  simulation,  the  results  should  show  the  degree 
of  degradation  in  filter  performance.  Finally,  to  compensate  for 
the  effects  of  dark  current  and  to  improve  tracking  performance, 
different  methods  in  changing  the  filter  structure  are  exploited. 

Mean  and  Variance  Convergence 

To  determine  the  appropriate  number  of  sample  runs,  two 
plots  are  generated.  One  plot  shows  the  sample  mean  error  versus 
the  number  of  runs,  and  the  other  plot  shows  the  sample  variance 
of  error  process  versus  the  number  of  runs.  These  plots  are  shown 
in  Fig.  (11)  for  the  Binomial  Method.  On  each  of  the  two  plots 
are  four  curves  corresponding  to  four  arbitrary  time  points 
(as  indicated  in  the  figure)  during  the  simulation  time. 

By  the  weak  law  of  large  numbers,  the  mean  and  variance  of  the 
error  process  should  converge  to  a steady  value  such  that  the 
variance  or  mean  of  the  error  process  does  not  change  much  in 
magnitude  with  each  addition  of  sample  runs.  Based  on  this 
reasoning  and  the  results  from  the  plots,  fifty  sample  runs  are 
selected.  This  number  of  runs  will  also  be  used  in  the  conditional 
method.  However,  it  does  not  appear  that  fifty  runs  are  sufficient 
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to  ensure  convergence  for  all  the  test  cases  listed  in  Tables 
1 and  2.  That  is,  the  performance  parameters,  X and  Y,  which  are 
defined  in  Chapter  II  as 


and 


X 


/2rr 

o 


(^TqJA  average  number  of  electrons 
occurring  in  one-half  the 
coherence  time  of  spot's 
dynamics  (A^=0)  (63) 


spot's  jitter  normalized  by  the  square  of 
the  beam  width  (64) 


f 


have  different  values.  The  plots  shown  in  Fig.  (12)  provide 
an  additional  illustration  of  confidence  in  the  sample  variances. 
They  are  plots  of  the  variances  of  the  error  process  (calculated 
from  sample  statistics)  for  different  numbers  of  sample  runs. 

On  the  true  variance  curves  are  symbols,  denoted  as  "x",  to 
distinguish  from  the  filter  variance  curve.  Along  with  the  curves 
are  the  upper  and  lower  bound  on  MSE  for  the  given  set  of  parameters 
X and  Y.  These  plots  show  that  as  the  number  of  simulation  runs 
increases,  the  fluctuations  of  the  true  and  filter  variances 
decrease.  The  plots  also  indicate  that  the  fluctuations  should 
settle  to  a value  that  is  between  the  upper  and  lower  bounds  (as 
expected).  Hence,  the  number  of  runs  indicate  the  amount  of 
confidence  that  the  sample  and  average  filter  sample  \ariances 
approach  the  corresponding  true  variances.  Although  fifty 
simulation  runs  may  not  be  appropriate  for  all  the  test  cases, 
it  does  provide  an  indication  of  filter  performance. 
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ERROR  VARIANCE 


Figure  12a  ^-POSITION  VpfflPINCES 
( 20  SAMPLE  RUNS) 
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ERROR  VARIANCE 
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Binomial  Method  (without  dark  current) 

For  the  Binomial  Method,  it  was  mentioned  in  the  previous 
chapter  that  simultaneous  events  can  be  generated  at  the  sample 
period.  At.  This  situation  would  occur  often  if  the  observation 
area,  AtAi.,  is  not  sufficiently  small.  . In  addition,  a low  rate 
function,  X(t,r),  is  selected  to  observe  the  propagation  of 
the  conditional  mean  and  variance  of  the  Snyder  Filter  and 
to  reduce  the  amount  of  computer  resources  by  this  method.  During 
this  part  of  the  study,  it  is  found  that  simultaneous  events  did 
occur  even  for  a rate  function  of  10  events  per  second  with  a 
sample  period  of  0.01  second.  For  X=5,  simultaneous  events 
occurred  on  the  average  of  two  to  three  times  per  sample  run. 

To  avoid  this  problem,  a software  routine  is  developed  to  randomly 
pick  one  of  two  simultaneous  events.  Note  that  this  routine 
ruins  the  statistical  properties  of  the  point  process. 

Figs.  (13a,  l4a)  and  (13b,  l4b)  are  plots  of  the  sample 
error  and  sample  filter  variance,  respectively,  for  the  parameters 
X=5  and  Y=l.  Note  on  the  time  axis  of  each  plot  are  the  event  times. 
In  Figs.  (13a)  and  (l4a),  the  errors  decrease  as  events  appear 
frequently  in  time  (as  expected).  Similarly,  in  Figs.  (13b)  and 
(l4b),  the  conditional  variance  is  updated  at  the  event  times, 
otherwise  the  variance  grows  in  time  when  an  event  is  not  observed. 
Note,  the  filter  variance  is  "driven"  to  a smaller  value  as  more 
events  appear  frequently  in  time.  If  the  filter  does  not  observe 
many  measurements,  the  filter  will  weigh  the  next  measurement 
heavily^ driving  the  filter  variance  to  a lower  value.  If  many 
previous  measurements  appear  frequently  in  time,  this  implies 
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a smaller  filter  growth  in  the  variance,  and  the  next  measurement 
will  have  a lower  weighting. 


For  this  research,  the  coherence  time,  t , of  the  spot 

is  one  second;  for  X=5,  the  term  a in  Eq.  (63)  is  equal 

n i o 

to  10  events  per  second.  This  is  the  average  rate  function  of 
the  Gaussian  spot,  after  integrating  X (t,r),  defined  in  Eq . (14), 

o 

over  the  detector  length  which  approximates  an  infinite  detector 
length  with  respect  to  the  Gaussian  beam.  It  is  well  known  for 
a poisson  process  that  the  mean  and  variance  of  the  counting 
process  N(t),  defined  in  Appendix  A,  are  respectively 
E[N(t)]  = /fT  X ( a , 6 ) dadf3  and  VAR[N(t ) ] = / fT  X(ct,e)dad3.  It 

AO  Vo 

then  follows  for  T=4  that  X (t,r)  can  be  expressed  as 


Ag(t,r)  = 


n I 


o/2tT 


/2tt 


exp( 


- (r--3  ( t ) ) ' 


2a 


10 


exp(zt£i|(iu:) 

/2ir  a 2a  d 


(64) 


so  that E[N  ( t )]  =40  events  with  standard  deviation  of  6.3  events. 
In  Figs.  (13)  and  (14),  the  number  of  events  for  this  particular 
case  are  34  and  36  events,  respectively.  Also,  in  these  sample 
runs,  simultaneous  events  did  not  occur.  To  the  second  moment 
level,  these  results  indicate  some  validity  in  the  Binomial 
Method  in  generating  the  point  process  observations  in  time  and 
space . 

Finally,  the  test  cases  for  this  method  are  plotted  in 
Figs.  (15)  through  (17);  note  the  different  scales  on  each  plot. 
As  expected.  Increasing  Y (the  spot's  jitter  normalized 
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by  the  square  of  the  beam  width)  increases  the  error  variances  - 
implying  poorer  performance.  Even  when  the  filter  is  initialized 
with  its  state  estimate  not  equal  to  the  true  position  of 
maximum  intensity  (i.e.,  initialize  from  non-zero  initial  conditions), 
the  plot  in  Fig.  ( 1 8 ) shows  that  the  error  process  is  zero-mean. 

This  is  expected  since  the  conditional  estimate  of  x(t)  is  unbiased 
and  x(t)  has  zero  mean. 


Conditi jnal  Method  (without  dark  current) 

Shown  in  Fig.  (19)  is  a sample  of  the  error  process  and 
the  filter  variance  generated  by  the  Snyder  Filter  for  the  set 
of  parameters,  X=20  and  Y=l.  From  these  plots,  it  is  evident 
that  the  measurements  of  event  locations  at  the  event  times 
are  important  in  minimizing  the  filter's  error  and  uncertainty. 
As  shown  in  Fig.  (19a),  the  errors  are  within  one  beam  width; 
and  in  Fig.  (19b),  there  is  little  filter  variance  growth 
as  events  appear  frequently  in  time  (as  seen  in  the  Binomial 
Method) . 


For  the  test  cases  in  Table  2,  it  is  desirable  to  choose 
values  of  the  sample  and  average  filter  variances  at  selected 
time  points.  For  each  case,  these  variance  points  are  plotted 
with  the  upper  and  lower  bound  curves.  However,  it  was  noted 
earlier  that  there  are  fluctuations  in  the  variances  which  make 
it  difficult  to  select  these  points.  A "time-average"  technique 
is  used  to  select  and  to  plot  these  points  easily.  This 
technique  is  depicted  in  Fig.  (20).  In  this  case,  a fixed  time 
interval,  that  is  chosen  to  be  much  greater  than  the  sample  period, 

contains  time  points  that  are  equally  spaced.  Corresponding 
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Figure  19a  X-P0S1T10N 
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Figure  19b  X-P0S1TI0N  FILTER  VRRlflNCE 
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Figure  20  "Time  - average"  Technique  (h(t)  denotes  time 

history  of  either  true  or  filter  variance  calculated 
from  sample  statistics) 


to  each  time  point  are  the  values  of  the  variances.  These  values 
are  then  summed,  and  the  resulting  sum  is  divided  by  the  number 
of  points  in  that  fixed  interval.  This  final  value  will  then 
be  used  to  denote  the  value  of  the  variance  at  the  midpoint  of 
the  fixed  time  interval.  This  time  interval  "runs"  across  the 
time  axis  to  calculate  these  "time-average"  values  until  the  end 
of  the  observation  time,  T. 

With  this  technique,  the  sample  and  average  filter  variances 

are  listed  in  Table  3 at  the  time  point,  t.  = 1.5  second.  These 

J 

values  are  plotted  with  the  upper  and  lower  bound  curves  in 
Figs.  (21).  The  time  point,  tj  = 1.5  second,  is  chosen  after 
examining  the  plots,  similar  to  plots  in  Fig.  (22),  for  each 
case  which  indicates  steady-state  is  reached.  Note  that  the  values 
of  the  true  and  filter  variances  are  close  to  the  upper  bounds 
which  seems  reasonable,  since  the  values  of  X=l,  X=5,  Y=l,  and 
Y=5  are  relatively  poor  conditions  for  good  tracking  performance. 
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Performance 

Parameter 

Normalized 

Bounds  on  MSE 

Normalized 
Simulation  Values 

X 

Y 

P* 

7? 

P* 

7? 

True 

Variance 

Filter 

Variant 

5.0 

0.5 

0.232 

1 

0.250 

0.233 

0.251 

5.0 

1.0 

0.358 

0.408 

0.388 

0.406 

5.0 

5.0 

0.905 

1.305 

1.340 

1.270 

5.0 

10.0 

1.320 

2.240 

2.370 

2.180 

1.0 

0.5 

0.366 

0.390 

0.340 

0.389 

1.0 

1.0 

0.618 

0.707 

0.621 

0.703 

1.0 

5.0 

1.790 

2.870 

2.570 

2.845 

1.0 

10.0 

2.900 

5.420 

4.900 

5.380 

0.5 

5.0 

2.320 

3.590 

4.200 

3.490 

10.0 

5.0 

0.659 

0.880 

0.935 

0.854 

15.0 

5.0 

0.5^5 

0.698 

0.670 

0.654 

20.0 

5.0 

0.496 

0.542 

0.592 

0.558 

0.5 

1.0 

0.732 

0.816 

0.898 

0.803 

10.0 

1.0 

0.270 

0.302 

0.319 

0.298 

15.0 

1.0 

0.227 

0.250 

0.243 

0.242 

20.0 

1.0 

0.200 

0.218 

0.222 

0.212 

Table  3.  Test  Cases  for  Conditional  Method 
(without  dark  current) 
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Figure  21b  Simulation  Values  of  True  and  Filter  Variances 
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Figure  22b  VRRiflNuE  CONVERGENCE 
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Figure  22c  X-P3SIT10N 
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That  is,  these  results  indicate  that  actual  performance  can 
be  described  by  the  upper  bound  on  MSE. 


Two  issues  should  be  of  concern  when  plotting  these  values. 
One  is  that  the  filter  performance  must  be  in  steady  state,  and 
the  other  is  that  the  number  of  sample  runs  in  the  simulation 
must  be  large  enough  to  yield  some  confidence  of  inferring 
information  from  the  error  variance  generated  by  either  sample 
statistics  or  averaging  of  filter  outputs.  Hence,  a particular 
case  is  chosen  to  extend  the  simulation  time  to  40  seconds;  the 
results  in  Fig.  (23)  indicate  that  steady  state  is  reached  in 
one  second.  Then  the  issue  again  is  of  confidence  in  the 
values  of  the  true  variances  which  can  be  increased  by  performing 
more  simulation  runs;  because  of  time  constraints,  fifty  simula- 
tions remain  unchanged.  To  illustrate  that  the  number  of  simula- 
tion runs  is  not  appropriate  in  all  cases,  refer  to  plots  in 
Fig.  (22).  For  this  particular  case,  X=20  and  Y=l,  the  mean  error 
is  not  converging  to  a definite  value. 

In  summary,  it  appears  that  the  upper  bound  describes  the 
actual  performance  of  the  Snyder  Filter  for  few  signal  counts 
per  coherence  time. 


Conditional  Method  (with  dark  current) 

Particular  cases  were  selected  to  examine  the  effects  of 
dark  current  and  to  determine  the  filter's  sensitivity  in 
performance  as  it  processes  these  noise  measurements.  The 
results  of  this  analysis  is  that,  for  average  signal-to-noise  counts 
ratio  of  500,  the  filter  was  very  sensitive  to  dark  distributed 
measurements.  The  results  are  plotted  in  Figs.  (24)  through 
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poorly  If  it  processes  even  one  dark  measurement.  Note  that 
the  plots  in  these  figures  are  not  close  to  ensemble  average 
statistics.  That  is,  to  obtain  the  sample  variance,  it  will 
require  a substantial  amount  of  additional  sample  runs.  Again, 
the  confidence  in  the  sample  statistics  is  an  issue.  On  the 
other  hand,  these  results  do  provide  insight  into  the  sensitivity 
of  the  Synder  Filter  to  dark  measurements.  However,  to  evaluate 
individual  effects  of  dark  current.  Figs.  (24a,  24b),  (25a,  25b), 
(26a,  26b),  and  (28a,  28b)  are  used  to  determine  which  sample 
runs  contain  dark  measurements.  Thus,  Figs.  (29a,  30a,  31a,  32a) 
and  (29b,  30b,  31b,.  32b)  are  plots  of  the  sample  errors  and  sample 
filter  variance,  respectively.  Figs.  (29)  and  (30)  are  for  average 
signal-to-noise  counts  ratio  of  500  and  50,  respectively,  for 
X=5  and  Y=l.  Similarly,  Figs.  (31)  and  (32)  are  for  average 
signal-to-noise  counts  ratio  of  500  and  5,  respectively,  for 
X=20  and  Y=l.  These  results  clearly  show  that  the  Snyder 
Filter  is  very  sensitive  to  dark  measurements.  Thus,  the 
Snyder  Filter  cannot  distinguish  between  a dark  or  signal 
measurement  (i.e.,  it  treats  all  events  as  signal). 

To  enhance  filter  performance,  filter-tuning  (in  this  case, 
increasing  the  beam-width)  and  residual  monitoring  were 
considered.  In  filter  tuning,  the  beam  width  is  increased, 
which  corresponds  to  weighing  the  measurements  less  heavily 
from  the  photodetector  surface.  In  other  words,  the  filter  is 
having  less  confidence  in  its  measurements.  Conversely,  if 
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Figure  26b  VRRlflNCE  CONVERGENCE 


L 


117 


ERROR 


119 


ERROR  + ONE  STANDARD  DEVIATION 


Figure  27d  X-POSITTON  VARIANCES 
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Figure  28c  X-P0SIT10N 


126 


127 


(WOlO'yj  N O 1 1 1 9 O d 


Figure  31a  X-POSITION 
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Figure  31b  X-POSHION  FILTER  VRRTpNCE 
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Figure  32b  X-POSITION  FILTER  VAR  I RNCE 


135 


the  beam  width  is  decreased,  which  means  that  the  signal  events 
are  distributed  at  locations  very  near  the  position  of  maximum 
intensity,  then  the  filter  should  weigh  the  measurement  heavily 
since  there  is  less  "corruption"  in  the  measurements.  Recall 
that  the  filter  assumes  the  photodetector  array  is  very  large; 
and  for  dark  limited  conditions,  the  dark  events  can  occur  any- 
where in  the  array  with  equal  probability.  Thus,  the  beam  width 
must  be  very  wide  in  order  to  put  less  weight  in  the  dark  measure- 
ments. Since  the  filter  cannot  distinguish  between  a dark  or 
signal  event,  the  signal  measurements  will  also  have  less  weight. 

In  other  words,  the  beam  is  assumed  to  be  wide  in  order  to  "wash" 
out  the  dark  events;  and  the  signal  events  are  also  assumed  to 
be  distributed  at  larger  distances  from  the  spot's  centroid  of 
maximum  intensity.  From  this  discussion,  it  becomes  evident  that 
altering  the  beam  width  in  order  to  compensate  for  the  dark 
events  has  some  problems,  especially  if  the  photodector  array  is 
dark  current-limited.  It  would  be  preferable  to  try  to  distinguish 
signal  events  from  dark  events  and  weigh  them  properly. 

In  residual  monitoring,  the  sequence  of  residuals  are 
monitored  in  real-time  to  decide  whether  the  current  measurement 
value  can  be  discarded  as  a "bad"  data  point.  Consider  Eq.  (23) 
which  is  repeated  below  as 

x(tk+)  = xCtk->  + K(tk-)[rk-x(tk-)]  (65) 

residual 

where  H(tk)=l.  From  Eq.  (65),  a residual  is  simply  the 

difference  between  the  measurement,  r,  , and  the  conditional 

* k5 


136 


A 

mean  xCt^.-)  of  the  measurement^  conditioned  on  only  the 
previous  history  of  measurements.  Residual  monitoring  would  be 
used  to  test  for  events  that  are  not  due  to  the  spot.  If  the 
events  fail  the  test,  then  the  filter  does  not  process  the 
dark  measurement.  Instead,  the  filter  propagates  its  conditional 
mean  and  covariance  until  it  senses  another  event  which  passes 
the  test.  For  this  study,  if  the  residual  in  Eq.  (23)  is 
greater  than  three  times  the  beam  width,  then  the  event  is 
discarded.  Note  that  in  residual  monitoring,  there  is  no  need 
to  increase  the  beam  width.  In  other  words,  the  event  locations 
(due  to  the  spot)  have  the  same  distribution  about  the  spot's 
centroid,  as  before.  Thus,  residual  monitoring  can  be 
exploited  effectively  for  high  signal-to-noise  counts  ratio, 
since  few  events  will  fail  the  test.  Physically,  residual 
monitoring  is  a better  method  of  compensation  than  filter 
tuning  in  order  to  handle  the  dark  current  effects. 

As  a rough  comparison  in  performance  improvement  between 
the  above  methods,  the  square  of  the  mean  errors  at  the  sample 
times  are  all  summed.  The  results  from  this  method  are  listed 
in  Tables  3 and  4.  The  tables  indicate  that  residual  monitoring 
is  a better  method  to  handle  the  dark  current  effects  than 
altering  the  beam  width.  For  particular  cases,  plots  of  the 
variances  in  which  these  methods  are  employed  are  shown  in  Figs. 
(33)  through  (36).  These  plots  show  that  the  Snyder  Filter 
performance  can  be  enhanced  substantially  through  residual 
monitoring  by  not  processing  the  dark  measurement.  As  expected, 
these  results  correspond  to  the  physical  arguments  mentioned  in 
previous  paragraphs. 
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Average 

Signal-to-Noise 
Counts  Ratio 


Method  of 
Compensation 


Beam  Width 
(xO.Ol  cm) 


Sum  of  the 
Square  of  the 
Mean  Error 


None 


5.0 


207.0 


500 

None 

5.0 

500 

Residual  Monitoring 

5.0 

500 

Filter  Tuning 

7.5 

500 

Filter  Tuning 

10.0 

500 

Filter  Tuning 

15.0 

500 

Filter  Tuning 

20.0 

500 

Filter  Tuning 

30.0 

50 

None 

5.0 

50 

Residual  Monitoring 

5.0 

50 

Filter  Tuning 

10.0 

50 

Filter  Tuning 

20.0 

50 

Filter  Tuning 

40.0 

50 

Filter  Tuning 

20.0 

1 

None 

5.0 

1 

Residual  Monitoring 

5.0 

1 

Filter  Tuning 

30.0 

1 

Filter  Tuning 

40.0 

555.0 

80.2 

323.0 

223.0 

166.0 

167.0 

196.0 

3150.0 

66.3 

1230.0 

198.0 

110.0 

377.0 

2860.0 

33.0 

350.0 
70.7 


Table  4.  Test  Cases  for  X=5  and  Y=1 
(with  dark  current) 
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Average 

Signal- to-Noise 
Counts  Ratio 

Method  of 
Compensation 

Beam  Width 
(xO.Ol  cm) 

Sum  of  the 
Square  of  the 
Mean  Error 

oo 

None 

5.0 

C'- 

C\J 

■=r 

500 

None 

5.0 

171.3 

500 

Residual  Monitoring 

5.0 

33.6 

500 

Filter  Tuning 

6.0 

142.4 

500 

Filter  Tuning 

7.0 

122.3 

500 

Filter  Tuning 

8.0 

95.1 

500 

Filter  Tuning 

10.0 

83.7 

500 

Filter  Tuning 

20.0 

57.8 

500 

Filter  Tuning 

30.0 

50.9 

50 

None 

5.0 

50 

Residual  Monitoring 

5.0 

41.0 

50 

Filter  Tuning 

10.0 

423.0 

50 

Filter  Tuning 

15.0 

208.0 

50 

Filter  Tuning 

30.0 

97.2 

50 

Filter  Tuning 

40.0 

105.0 

5 

None 

5.0 

8940.0 

5 

Residual  Monitoring 

5.0 

60.6 

5 

Filter  Tuning 

20.0 

882.0 

5 

Filter  Tuning 

30.0 

396.0 

Table  5.  Test  Cases  for  X=20  and  Y=1 
(with  dark  current) 
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Figure  33a  /-POSITION  VaRlANCES 
(RESIDUAL  MONITORING) 
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Figure  35b  X-P0SIT1ON  '/HR I HNCJTS 
(FILTER  TUNING) 


ERR3R  VnRI PNCE 


Figure  36a  X-P0SIT10N  VRR1RNCES 
(RESIDUAL  MONITORING) 


Figure  36b  X-POSITION  VRRIRNCES 


(FILTER  TUNING) 
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Summary 


For  the  test  cases  considered,  the  filter's  average 
performance  lies  near  the  upper  bound  curve,  which  seems  reasonable 
since  the  set  of  parameters  is  selected  to  yield  relatively 
poor  tracking  (i.e.,  few  signal  counts  per  coherence  time). 

When  dark  current  is  included  in  the  measurement  model 
of  the  Snyder  Filter,  the  estimator  is  found  to  be  very 
sensitive  to  dark  measurements  and  tracking  performance  is 
degraded  substantially. 

Based  on  physical  arguments  and  simulation  results,  residual 
monitoring  provides  a better  method  to  compensate  for  the  effects 
of  dark  current  than  filter  tuning. 
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V Conclusions  and  Recommendations 


i, 


Conclus ions 

Assuming  no  dark  current,  the  average  error  performance 
of  the  Snyder  Filter  can  be  described  by  steady-state  upper 
and  lower  bounds  on  mean  square  error  (MSE).  These  bounds  are 
evaluated  as  a function  of  two  physically  motivated  parameters. 
These  parameters  are  described  as:  (1)  the  average  number  of 
signal  photoelectrons  observed  in  a coherence  time  of  spot's 
dynamics,  and  (2)  the  spot's  jitter  normalized  by  the  square 
of  the  beam  width.  Based  on  these  parameters,  sufficient 
conditions  are  presented  such  that  the  bounds  converge  to  the 
actual  performance.  When  sufficient  conditions  are  not  estab- 
lished, the  results  from  the  Monte  Carlo  simulation  indicate 
that  actual  performance  can  be  described  by  the  upper  bound 
(in  this  study,  the  conditions  correspond  to  low  signal  counts 
per  coherence  time). 

Presenting  the  dark  measurements  in  the  simulation  gave 
further  insights  into  the  tracking  capabilities  of  the  Snyder 
Filter.  Even  for  an  average  signal-to-noise  counts  ratio  of 
500,  the  results  show  that  this  estimator  was  highly  sensitive 
to  dark  measurements.  Thus,  the  Snyder  Filter  treats  all  events 
as  signal. 

Based  on  simulation  results  and  physical  arguments,  residual 
monitoring  is  a better  way  to  handle  the  dark  events  rather  than 
increasing  the  beam  width.  This  can  be  explained  as  follows: 
by  increasing  the  assumed  beam  width,  the  Snyder  Filter  puts  less 
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weight  on  the  dark  measurements.  At  the  same  time,  more  events 
will  be  assumed  to  be  distributed  from  the  spot's  centroid. 

Thus,  the  filter  will  not  weigh  heavily  the  measurement  due  to 
the  spot's  Gaussian  intensity  profile.  Also,  the  Snyder  Filter 
assumes  the  measurements  are  from  a large  photosensitive  surface, 
and  therefore  the  beam  width  must  be  very  wide  to  "wash"  out  the 
dark  current  effects.  On  the  other  hand,  residual  monitoring  is 
a method  that  involves  testing  for  dark  events.  In  this  study, 
if  a measurement  (i.e.,  event  location)  occurred  that  is  three 
beam  widths  away  from  the  position  of  maximum  intensity,  the  Snyder 
Filter  would  not  process  this  measurement.  Hence,  it  is  unnecessary 
to  increase  the  beam  width  in  residual  monitoring.  For  high  signal- 
to-noise  counts  ratio,  residual  monitoring  can  be  used  effectively. 
Results  from  the  simulation  demonstrate  this  physical  reasoning. 

Recommendations 

An  area  of  recommended  research  is  to  conduct  this  study  in 
two  dimensions  which  would  investigate  different  equal  intensity 
ellipsoidal  shapes  of  the  spot's  intensity  profile.  Also,  investi- 
gation of  other  intensity  profiles  besides  the  assumed  Gaussian- 
shaped  may  warrant  further  research.  The  effects  in  filter 
performance  can  then  be  investigated. 

Another  issue  of  research  is  how  to  incorporate  the  point- 
process  measurements.  This  study  was  performed  in  which  the 
rate  function  was  chosen  to  have  a low  value,  but  in  the  real 
world,  rates  of  10^  events  per  second  can  occur  on  the  detector 
surface.  This  imposes  a very  strict  requirement  of  computation 
on-line.  A suggestion  was  made  to  store  the  photoelectrons  in 
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charged-coupled  devices.  Then  the  Snyder  Filter  can  process 
its  measurement  at  fixed  time  intervals.  In  this  case,  multiple 
measurements  are  available  to  the  Snyder  Filter.  Then,  how 
does  one  process  these  stored  or  "batch"  measurements,  and  how 
does  one  account  for  different  arrival  times  (if  at  all)?  This 
creates  an  area  of  research  that  warrants  further  investigation. 

Another  illustration  in  how  to  Incorporate  the  point-process 
measurements  is  the  issue  of  simultaneous  events  encountered  in 
the  study.  How  would  one  process  these  multiple  measurements 
at  single  time  instant  in  the  structure?  Obviously  there  is  more 
information  concerning  the  spot's  position  of  maximum  intensity 
(assuming  no  dark  current).  This  question  should  be  answered 
with  further  research. 

This  study  analyzed  the  Snyder  Filter  in  open-loop  configura- 
tion; thus,  a follow-on  study  with  control  inputs  to  position- 
sensitive  devices,  such  as  a mirror,  warrants  further  investiga- 
tion . 

A final  issue  and  recommendation  is  to  develop  ad-hoc 
methods  in  how  to  handle  the  dark  current  effects  for  low  signal- 
to-noise  counts  ratio.  Since  residual  monitoring  was  found  to 
be  a more  effective  method  of  compensation  than  filter  tuning, 
different  methods  using  residual  monitoring  should  be  investigated. 
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Appendix  A 


Statistics  of  a Time-Space  Point  Process 

The  purpose  of  this  appendix  is  to  derive  the  statistics 
of  a point-process  in  time  and  space  based  on  an  incremental 
description.  In  its  application  of  this  report,  we  define 
a rate  function  to  be 

X(t ,r)  = JA  |Ud(t,r)  1 2dr;  dr  = dxdy 

1 o d 

where  t denotes  the  time  continuum,  r is  the  dummy  variable 
representing  a location  on  a photodetector  surface,  q is  the 
quantum  efficiency,  h is  Planck's  constant,  A^  is  the 
detector  surface  area,  is  the  complex  field  envelope,  fQ 
is  frequency  of  the  input  field,  dx  is  differential  length 
along  x-axis,  dy  is  differential  length  along  y-axis  and  |*| 
is  modulus  of  field  envelop.  This  rate  function  describes  the 
intensity  of  events  appearing  on  the  detector  surface  in  time 
and  space. 

Shown  in  Pig. A1  is  a detector  array  of  photon  or  quantum 
detectors  which  can  measure  the  temporal  and  spatial  location 
of  a photon-to-electron  conversion. 

It  is  assumed  that  each  photon-to-electron  conversion  or 
event  on  the  detector  surface  is  independent.  Then  based  on 
its  Incremental  properties  the  x-axis,  y-axis,  and  t-axis  are 
partitioned  into  intervals  as  shown  in  Fig.  A-2: 

Then,  with  the  incremental  description  is  as  follows 
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p 


event  occurring 
in  AtAA  at 


X(t,r)AtAA  + o(AtAA)  (A-l) 


P 


no  event  occurring 
in  AtAA  at 


1 -A ( t , r ) AtAA  + o(AtAA)  (A-2) 


P 


more  than  one 
event  in  AtAA 


<W 


o ( At AA ) 


( A — 3 ) 


vjhere  P[*]  denotes  probab ility,  of  o(AtAA)  are  higher  order 

terms  of  o(AtAA)  with  limit  = 0.  It  is  desirable  to 

AtAA+0 

obtain  a statistical  description  of  event  times  t^'s  event 

locations  r^'s,  and  the  number  of  events  n,  appearing  on  entire 

detector  surface  A^,  in  a time  interval  [0,t].  Thus,  we  need 

the  joint  probability  density  function  of  these  random  variables 

t,  , r,  , and  n (i.e.  f[{t.  },{r,  },n]  where  {•}  denotes  set  of  event 
k k k k 

locations  or  times.  f[{t. },{r  },n]  completely  characterizes 

K rC 

the  output  process  from  the  detector  surface. 

To  obtain  f[  { t ^ } , {r  } ,n] , consider  a counting  process 
which  is  depicted  in  Fig.  A-3. 

The  point  process  probability  distribution  for  N^Ci.e., 

P[N,=  n])  will  be  derived  using  characteristic  functions.  Note 
that  the  process  is  the  number  of  events  regardless  of  when 
and  where  the  event  occurred.  Since  is  a discrete  random 
variable , 


jvN, 


E[e  = eJ<'0)  P[Nt  = 0]  + eJVP[>t  = 1] 

+ eJ’2vP[Nt  = 2]+  ...  + eJV  P[Nfc  = n] 


JVT 


15't 


(A-4) 


Fig  A. 3 Time-Space  Counting  Process 


1 rc. 


Let  be  a random  variable  that  is  equal  to  a one  in  which 

an  event  occurred  in  AtAA  or  a zero  in  which  case  no  event 

n 

occurred.  It  then  follows  that  N = £ n.  where  k = 0,1,2 ,3. 

k=l 

Then  Eq.  (A--'0  becomes 


E[e°  = E[exp( jv£n  ) ] 


= U E [exp ( j vn  ) ] 

k=l  K 


(A-5) 


where  we  have  use  the  property  of  statistical  independence  of 
events.  Upon  using  Eqs.  (A-l)  through  (A-3),  Eq.  (A-5)  becomes 


jvk.  Nt 


] = if  eJvlA(t.  , r,  ) AtAA  + ejv°[l-A(t  ,r  )AtAA]+  o ( At AA ) 

k=l  k K 


n A ( t , r ) AtAA(eJV-l)  +1+0  (AtAA) 
k=l  k K 


exp  l ln[  1 + A(t  ,r  )AtAA(eJV-l)  + o(AtAt)]  (A-6) 

k=l  R K 


and  using  ln(l+z)  ; z for  small  z,  Eq  (A-6)  becomes 
jvN  Nt 

E[e  fc]  = exp  [ l A ( t ,r  ) At AA ( V-1 )] 
k=l 


Also,  as  At  A A->-0 , the  sum  in  the  above  relation  becomes  a definite 


integral  so  that 

j vN  . t — — 

E[e  t]  = exp[(eJV-l)  J T X(a,a)d«dB] 


( A-7 ) 


A . 0 
a 


with  the  series  expansion  of  exp(x)  = 1 + x + + 
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Eq  (A-7)  yields 


jvN 

E[e  ]=  exp  ( -vi . ). 

a ,t 


e2j' v[exp(p  )]2 

d,t 


wnere  p 


+ . . . + e^v[exp(p  )] 
Hd  ,t 

n! 


= / ) X(t,r)da,d3 


) 


(A-8) 


djt  A 0 
a 


Comparing  Eq.  (A-8)  with  Eq.  (A-4),  we  obtain 


P[Nt  = n]  = — ^ 


(p  )nexp  (-p  ) 

a ,t 


n! 


(A-9) 


Eq  (A-9)  describes  the  counting  process  in  interval 
[0,t]  without  regard  to  the  event  times  and  event  locations. 

To  obtain  f[{t^},  (rk},n],  and  using  independence  of  events, 
consider  the  following  probabilities: 


P[one  event  at  (t^i^)]  = X(t1,r1)AtAA  (A-10) 

P[one  event  at  (t1,r1)  and  no  events  in  [t,T]]= 

= X(t  ,r  )AtAA  exp(-  / /t  X(a,‘3)dad3  (A-ll) 

A 4- 
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Finally,  the  probability  density  function  of  the  ordered 
event  times,  event  locations,  and  the  number  of  events  in 
[0,T]  on  A^  is  defined 


f [{tkJ , {rk} ,n] 


lim 


At  AA-*0 


pLd 

( At AA)n 


( A-12 ) 


where  P[*]  denotes  the  probability  of  one  event  at  (t^,r^),  one 
at  (t^r^),  - . . , one  at  (t  ,r^).  Using  the  Eqs  (A-10)  and 
(A-ll),  (A-12)  will  yield 

i-  - - n 

f[(t,  },{r  },n]  = exp  ( - / / A ( a , g )dad g) n X ( t ,r,  ) (A-13) 

k k A,  0 k=l  k K 

d 


Eq  (13)  completely  characterizes  the  statistical  output  process 
of  the  detector  surface. 

Since  this  report  simulates  the  detector  output  in  one 
dimension,  simply  replace  r^  with  r^,  to  yield 

t n 

f{tkJrk,n}  = exp(“  / I / A(a,g)(a,g)  n A(tk,tk)  (A-l4) 

0 k — 1 

where  i is  detector  length. 

To  obtain  the  joint  probability  density  of  event  times  and 
number  of  events  without  regard  to  where  the  events  occurred 
(i.e.  f{{tk),n}),  Eq  (A-ll)  will  yield  the  marginal  density 

t n 

f({t.  },n)  = exp(-  / / A(oi,g)(a,g)  n / A(a,g)da 

K x.  o i=l  * 

where  the  rate  function  % ^rkstk^  integrated  over  the 

detector  length  yielding  no  spatial  dependence  on  j A(a.  ,t )da 

*,  k k 
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This  completes  the  statistical  description  of  a time-space 


point  process  which  characterizes  the  output  of  a photo- 
sensitive surface. 


Appendix  B 


Upper  and  Lower  Bounds  on  Mean  Square  Error 

i1 

The  Snyder  Filter  is  a minimum  mean-square-error  estimator 
that  uses  time-space  point  process  observations.  The  purpose 
of  this  appendix  is  to  solve  the  upper  and  lower  bounds  on 
mean-square-error  (MSE)  in  order  to  evaluate  the  estimation 

performance  of  the  Snyder  Filter.  Also,  the  bounds,  derived 

* 

in  Ref.  4 are  solved  in  steady-state  (i.e.,  ^ = 0). 

dt  dt 

Parameters  are  defined  from  the  dynamical  equation  of  P*  and  P* 
to  provide  meaningful  physical  descriptions  about  the  actual 
performance  of  the  Snyder  Filter.  Based  on  these  parameters, 
sufficient  conditions  can  be  established,  such  that  the  bounds 
are  identical;  thus,  actual  performance  can  be  determined. 

However,  when  the  parameters  do  not  establish  sufficient 
condit ions , the  bounds  are  no  longer  identical,  and  actual 
performance  is  difficult  to  obtain.  These  parameters  will 
now  be  derived  in  one  dimension  which  involves  only  estimation 
performance . 

The  parameters  are  derived  in  one-dimension  for  a spot's 
dynamics  to  be  modelled  as  a Scalar  First-Order  Markov  process. 

In  addition,  this  solution  only  involves  estimation  performance. 

Consider  the  spot's  centroid  (or  position  of  maximum  intensity) 
which  is  modelled  by  a linear  stochastic  differential  equation 

dx( t ) = F(t)x(t)dt  + G ( t )d$( t ) ( B-l ) 

where  x(t)4  position  of  maximum  intensity,  6 ( t ) is  a Wiener 
process,  F(t)  is  the  state  coefficient,  and  G(t)  is  the  noise 

16Q 

A 


input  coefficient.  The  time-space  point  process  has  a rate 
function  modelled  as 


where 


X ( r , x ( t ) , t ) = u(t)Y(r,x(t)) 


nl  (t) 

p(t)  = —2 J?iTo 

hf 


2<j2 


(B-2) 


( B—  3 ) 


(B-'J) 


with  r is  a dummy  variable  for  spatial  location  on  detector, 
ar.d  a is  the  beam  width. 

The  lower  bound  on  MSE, denoted  as  P*(t),  is  described  by 
the  following  differential  equation  that  is  specialized  in 
one-dimensional  space  as 

dP# 

- F(t)Px(t)  + P#(t)F(t)  + G(t)Q(t)G(t) 


(B-5) 


where  u - E[u(t)];  E[«]  = expected  value  operator 
Q(t)  is  strength  of  6ft) 

The  lower  bound  is  specialized  for  the  following  assumptions: 

(1)  x is  modelled  as  an  output  process  of  a First  Order 

Lag;  then  F(t)  = -i  and  its  state  transition  function  is 

T c 

<*>  ( t , t ) = e ^ t)/tc.  The  strength  Q of  the  white  noise 

driving  the  First  Order  lag  is  given  in  Ref  (9)  as  Q =— — 0 

Tc  x 

161 


to  obtain  stationary  statistics  at  the  output.  a ' 
is  the  mean  square  value  of  the  spot's  displacement  on 

o -( t ) T 

the  detector  surface  (i.e.  E[x(t )x(t+x ) ] = a e 
and  (2)  G(t)  = 1.  Eq  (R-6)  will  yield. 


* - “2p  , n V p2 

..  W ' ~2  P* 
d t c a 


( B-6 ) 


Solving  for  steady  state  _ * 

dt 

equation  which  is  expressed  as 


= 0,  yields  a quadratic 


-221  + 
2tcp 


i [<  ♦ ,(9jL,]h 


■ (4  )(  • [ 

■ H X-  * [ 


1 * 

p 2a 


1 + !»(--%  ) (^£_) 


P T ~~  Ql 

Dividing  Eq  ( B— 8 ) by  a and  defining  X=  y = — % 


yields 


(B-7 ) 


■H 


(B-8 ) 


P*  J_  / 

7 " 2X  \ " 


1 + [1  + 4xy] 


!V]'S  ) 


(B-9 ) 


But , u 


/ 2 tt  a and  Q = — — 


, and  substituting 


into  x and  Y gives 


* - (hrr  To  < TT1 


Y = 


(B-10) 
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X describes  the  average  number  of  photons  detected  in  a coherence 
time  of  the  spot's  dynamics  and  Y results  in  the  spot's  jitter 
normalised  by  the  square  of  the  beam  width.  These  parameters 
are  discussed  in  detail  in  Chapter  II. 

The  relation  for  the  upper  bound  on  MCE  is  expressed  as 


* 

dP 

dt 


# O 1 

FP  + P*F  + GQG  - pP*[P*  + o ] P* 


(B-ll) 


Applying  the  same  assumptions  as  in  P„ , Eq  (B-ll)  yields 


dP!  = =1  P*  + o - 

dt  t J * 2 

c p +o 


dP* 

For  = 0 yields  another  quadratic  equation  which  solves  for 

dt 


1 

2 


1 

2 


^ CL 

dO 

Q - tc 

r 2a  2 1 

Q-  Tc  4c  0 Tp 

-L 

'2 

c 

2 , - 2 ± — 

T + P T + " 

C C J 

h 

Q,0-2<,2 

Qt  -2o^  MQo't 

C 1 c 

2+tttc 

2 + liTc  ?+^Tc 

( 1 + y t ) 
c 


/ Q.T  \ 

(!+”) 


&-) 


+ 4 


(£)(>♦*) 


(B-12) 


Again,  letting  X = 


UT 


Qt 


, 2 

^ , and  dividing  o 


P n 


into  Eq  (B-12)  yields 


* 

P 

T 

a 


(5T?fTT)((Y-1,?  + 


+ Mx+Dy)]’5 


where  X and  Y is  defined  in  (B-10).  Note  that  for  X >> 
Y<<1,  the  bounds  converge  to  actual  mean  square-error, 
conditions  are  discussed  in  Chapter  TI . 


1 
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